import numpy as np
from numpy import *
import scipy
from scipy.optimize import minimize 
from scipy.integrate import nquad
import matplotlib.pyplot as plt
import matplotlib
import time
import sys
from numba import jit
from sympy import re, im, I, E
from sympy import exp, N, S
from sympy.matrices import Matrix
matplotlib.use('Agg')

################################################################################################################################
#Use input 96 for agreement with exp work arXiv:2109.12631 (corresponds to D_0=0)
################################################################################################################################
# Model parameters - Using numbers in https://arxiv.org/pdf/1901.10485.pdf

#lattice spacing (in m, we take graphene *primitive* lattice constand to be order and angstrom)
a=1.42*10**(-10) 
# Dirac momentum of 1 layer of graphene
kD=4*np.pi/(3*a*sqrt(3))
# Twist angle
theta=1.53*np.pi/180
# momentum scale of MBZ - ktheta is length of one side of hexagon of MBZ
ktheta=2*kD*np.sin(np.abs(theta)/2)
# alpha=w_{AB}/(vF*kD*sin(theta)) - dimensionaless constant to characterize magic angle
alpha=.414
# fermi velocity of dirac cone in MBZ
vF=1/(kD*theta)
hbar = 6.6*10**(-13)
vF=10**6*hbar
# Note energies are of order vF*kD*theta

# kappa=wAA/WAB - assume wAB12=wAB23
phaseStep=12

kappa=int(int(sys.argv[1])/(phaseStep))*1/(phaseStep-1)
#inter-subllatice couplings - in mirror symmetric model they are equal for 12, 23
wab12=alpha*vF*kD*theta
wab23=alpha*vF*kD*theta
#intra-subllatice couplings - in mirror symmetric model they are equal for 12, 23 - determined by ratio kappa
waa12=kappa*wab12
waa23=kappa*wab23

#displacement vector. Set it to zero to achieve mirror symmetric limit
d=0
# strength of displacement field. We choose it to be of order of interlattice hopping

D=int(int(sys.argv[1])%phaseStep)/(phaseStep-1)*wab12*.5

print('w0 ',kappa,'interaction step: ',D/wab12)

################################################################################################################################
# Construction of MBZ. Define primitive lattice vectors in MBZ

q1x=0
q1y=-ktheta

q2x=ktheta*np.sqrt(3)/2
q2y=ktheta*1/2

q3x=-ktheta*np.sqrt(3)/2
q3y=ktheta*1/2

# Bravais lattice vectors in MBZ
b1x=np.sqrt(3)/2*ktheta
b1y=3/2*ktheta
b2x=np.sqrt(3)/2*ktheta
b2y=-3/2*ktheta

# Number of shells used to construct MBZ - ie add another "shell" of Moire unit cells nshells*b_reciprocal from gamma point of zeroth MBZ
nshells=3
radius=np.sqrt((nshells*b2x+q1x)**2+(nshells*b2y+q1y)**2)+.001*ktheta

shell1momentax=[]
shell1momentay=[]

shell2momentax=[]
shell2momentay=[]

shell3momentax=[]
shell3momentay=[]

shell4momentax=[]
shell4momentay=[]

shell5momentax=[]
shell5momentay=[]

for a in range(-10,10):
	for b in range(-10,10):
		if np.sqrt((a*b1x+b*b2x)**2+(a*b1y+b*b2y)**2)<np.sqrt((1*b2x)**2+(1*b2y)**2)+.001*ktheta:
			if a!=0 or b!=0:
				
				shell1momentax.append(a*b1x+b*b2x)
				shell1momentay.append(a*b1y+b*b2y)
for a in range(-10,10):
	for b in range(-10,10):
		if np.sqrt((a*b1x+b*b2x)**2+(a*b1y+b*b2y)**2)<np.sqrt((2*b2x)**2+(2*b2y)**2)+.001*ktheta and np.sqrt((a*b1x+b*b2x)**2+(a*b1y+b*b2y)**2)>np.sqrt((1*b2x)**2+(1*b2y)**2)+.001*ktheta:
			

			shell2momentax.append(a*b1x+b*b2x)
			shell2momentay.append(a*b1y+b*b2y)

for a in range(-10,10):
	for b in range(-10,10):
		if np.sqrt((a*b1x+b*b2x)**2+(a*b1y+b*b2y)**2)<np.sqrt((3*b2x)**2+(3*b2y)**2)+.01*ktheta and np.sqrt((a*b1x+b*b2x)**2+(a*b1y+b*b2y)**2)>np.sqrt((2*b2x)**2+(2*b2y)**2)+.001*ktheta:
			
			shell3momentax.append(a*b1x+b*b2x)
			shell3momentay.append(a*b1y+b*b2y)

for a in range(-10,10):
	for b in range(-10,10):
		if np.sqrt((a*b1x+b*b2x)**2+(a*b1y+b*b2y)**2)<np.sqrt((4*b2x)**2+(4*b2y)**2)+.001*ktheta and np.sqrt((a*b1x+b*b2x)**2+(a*b1y+b*b2y)**2)>np.sqrt((3*b2x)**2+(3*b2y)**2)+.001*ktheta:
			shell4momentax.append(a*b1x+b*b2x)
			shell4momentay.append(a*b1y+b*b2y)
for a in range(-10,10):
	for b in range(-10,10):
		if np.sqrt((a*b1x+b*b2x)**2+(a*b1y+b*b2y)**2)<np.sqrt((5*b2x)**2+(5*b2y)**2)+.001*ktheta and np.sqrt((a*b1x+b*b2x)**2+(a*b1y+b*b2y)**2)>np.sqrt((4*b2x)**2+(4*b2y)**2)+.001*ktheta:
			shell5momentax.append(a*b1x+b*b2x)
			shell5momentay.append(a*b1y+b*b2y)	
		

################################################################################################################################
# Moire Potential construction
phi=2*np.pi/3

Id2=np.block([[1,0],[0,1]])
zero2=np.block([[0,0],[0,0]])
paulix=np.block([[0,1],[1,0]])
pauliy=np.block([[0,-1j],[1j,0]])
pauliz=np.block([[1,0],[0,-1]])



T12_1=np.conjugate(np.transpose(np.block([[waa12,wab12],[wab12,waa12]])))*np.exp(1j*d/2*ktheta)
T12_2=np.conjugate(np.transpose(np.block([[waa12,wab12*np.exp(1j*phi)],[wab12*np.exp(-1j*phi),waa12]])))*np.exp(-1j*d/4*ktheta)
T12_3=np.conjugate(np.transpose(np.block([[waa12,wab12*np.exp(-1j*phi)],[wab12*np.exp(1j*phi),waa12]])))*np.exp(1j*d/4*ktheta)

T12_1c=np.conjugate(np.transpose(np.block([[waa12,wab12],[wab12,waa12]])))*np.exp(-1j*d/2*ktheta)
T12_2c=np.conjugate(np.transpose(np.block([[waa12,wab12*np.exp(1j*phi)],[wab12*np.exp(-1j*phi),waa12]])))*np.exp(1j*d/4*ktheta)
T12_3c=np.conjugate(np.transpose(np.block([[waa12,wab12*np.exp(-1j*phi)],[wab12*np.exp(1j*phi),waa12]])))*np.exp(-1j*d/4*ktheta)


T23_1=np.conjugate(np.transpose(np.block([[waa23,wab23],[wab23,waa23]])))*np.exp(-1j*d/2*ktheta)
T23_2=np.conjugate(np.transpose(np.block([[waa23,wab23*np.exp(1j*phi)],[wab23*np.exp(-1j*phi),waa23]])))*np.exp(1j*d/4*ktheta)
T23_3=np.conjugate(np.transpose(np.block([[waa23,wab23*np.exp(-1j*phi)],[wab23*np.exp(1j*phi),waa23]])))*np.exp(-1j*d/4*ktheta)

T23_1c=np.conjugate(np.transpose(np.block([[waa23,wab23],[wab23,waa23]])))*np.exp(1j*d/2*ktheta)
T23_2c=np.conjugate(np.transpose(np.block([[waa23,wab23*np.exp(1j*phi)],[wab23*np.exp(-1j*phi),waa23]])))*np.exp(-1j*d/4*ktheta)
T23_3c=np.conjugate(np.transpose(np.block([[waa23,wab23*np.exp(-1j*phi)],[wab23*np.exp(1j*phi),waa23]])))*np.exp(1j*d/4*ktheta)


# number of points on first diagonal of grid over first MBZ. nptsd should be an even number as there may be complications when the M point is included (none of our results have used an odd grid)
nptsd=10
npts=3*nptsd**2-3*nptsd+1-2*nptsd-(nptsd-2)
print('npts: ',npts)
domainsize = (3*ktheta**2*np.sqrt(3)/2/3)

# parameters used to define grid in k-space. scale can shrink the grid if you want to move the grid off of the high symmetry points. However, to construct a grid consistent with scattering, scale should be set to 1 (must be 1 for all HF calculations)
scale=1
vx=scale*ktheta*1/2
vy=scale*ktheta*np.sqrt(3)/2
c1x=scale*1/2*ktheta*1/(nptsd-1)
c1y=-scale*np.sqrt(3)/2*ktheta*1/(nptsd-1)
c2x=-scale*1/2*ktheta*1/(nptsd-1)
c2y=-scale*np.sqrt(3)/2*ktheta*1/(nptsd-1)
by=scale*0
bx=-scale*ktheta*1/(nptsd-1)
count=0
rq1x=3*ktheta/(2*np.sqrt(3))
rq1y=np.sqrt(3)/2*ktheta


# These are to store the momenta values of the G grid. 
qxvalstot=[]
qyvalstot=[]
# store layer index for each q point (layer 1 and 3 will share some)
qvalslayer=[]

for i in range(-10,10):
	for j in range(-10,10):
		if np.sqrt((i*b1x+j*b2x+q1x)**2+(i*b1y+j*b2y+q1y)**2)<radius:

			qxvalstot.append(i*b1x+j*b2x+q1x)
			qyvalstot.append(i*b1y+j*b2y+q1y)

			qvalslayer.append(2)
for i in range(-10,10):
	for j in range(-10,10):
		if np.sqrt((i*b1x+j*b2x+q1x+q2x)**2+(i*b1y+j*b2y+q1y+q2y)**2)<radius:

			qxvalstot.append(i*b1x+j*b2x+q1x+q2x)
			qyvalstot.append(i*b1y+j*b2y+q1y+q2y)

			qvalslayer.append(1)

for i in range(-10,10):
	for j in range(-10,10):
		if np.sqrt((i*b1x+j*b2x+q1x+q2x)**2+(i*b1y+j*b2y+q1y+q2y)**2)<radius:

			qxvalstot.append(i*b1x+j*b2x+q1x+q2x)
			qyvalstot.append(i*b1y+j*b2y+q1y+q2y)

			qvalslayer.append(3)
qxvalstotminus=[]
qyvalstotminus=[]
# store layer index for each q point (layer 1 and 3 will share some)
qvalslayerminus=[]

for i in range(np.size(qxvalstot)):
	qxvalstotminus.append(-qxvalstot[i])
	qyvalstotminus.append(-qyvalstot[i])
	qvalslayerminus.append(qvalslayer[i])

n=np.size(qxvalstot)*2

# rotate momenta counterclockwise by pi
def C2zrotate(kx,ky):
    kxprime=-kx
    kyprime=-ky
    return kxprime,kyprime
# rotate momenta counterclockwise by 2pi/3
def C3rotate(kx,ky):
    kxprime=np.cos(2*np.pi/3)*kx-np.sin(2*np.pi/3)*ky
    kyprime=np.sin(2*np.pi/3)*kx+np.cos(2*np.pi/3)*ky
    return kxprime,kyprime
# C2z symmetry operator in one spin flavor
def C2zfull():
	nbands=np.size(qxvalstot)*2
	UC2z=np.full((nbands*2,nbands*2),0.+0.j)
	for i in range(np.size(qxvalstot)):
		kxprime,kyprime=C2zrotate(qxvalstot[i],qyvalstot[i])


		x=np.where(np.round(qxvalstotminus/ktheta,3)==round(kxprime/ktheta, 3))
		y=np.where(np.round(qyvalstotminus/ktheta,3)==round(kyprime/ktheta, 3))

		c2zindex=0

		intersect=np.intersect1d(x, y)

		for j in range( np.size(intersect)):
			if qvalslayer[i]==qvalslayerminus[np.intersect1d(x, y)[j]]:
				c2zindex=np.intersect1d(x, y)[j]


		UC2z[2*i][nbands+2*c2zindex+1]=1
		UC2z[2*i+1][nbands+2*c2zindex]=1
		UC2z[nbands+2*c2zindex][2*i+1]=1
		UC2z[nbands+2*c2zindex+1][2*i]=1




	blockzeros=np.full((nbands,nbands),0)
	return UC2z
# C2z symmetry operator with Spin
def C2zSpin():
	nbands=np.size(qxvalstot)*2
	UC2z=np.full((nbands*2,nbands*2),0.+0.j)
	for i in range(np.size(qxvalstot)):
		kxprime,kyprime=C2zrotate(qxvalstot[i],qyvalstot[i])


		x=np.where(np.round(qxvalstotminus/ktheta,3)==round(kxprime/ktheta, 3))
		y=np.where(np.round(qyvalstotminus/ktheta,3)==round(kyprime/ktheta, 3))

		c2zindex=0

		intersect=np.intersect1d(x, y)

		for j in range( np.size(intersect)):
			if qvalslayer[i]==qvalslayerminus[np.intersect1d(x, y)[j]]:
				c2zindex=np.intersect1d(x, y)[j]


		UC2z[2*i][nbands+2*c2zindex+1]=1
		UC2z[2*i+1][nbands+2*c2zindex]=1
		UC2z[nbands+2*c2zindex][2*i+1]=1
		UC2z[nbands+2*c2zindex+1][2*i]=1





	blockzeros=np.full((2*nbands,2*nbands),0)
	UC2zFull=np.block([[UC2z,blockzeros],[blockzeros,UC2z]])
	return UC2zFull
# C2zT operator in one spin flavor and a single valley
def C2zT():
	nbands=np.size(qxvalstot)*2
	UC2z=np.full((nbands*2,nbands*2),0.+0.j)
	for i in range(np.size(qxvalstot)):


		UC2z[2*i][2*i+1]=1
		UC2z[2*i+1][2*i]=1
		UC2z[nbands+2*i][nbands+2*i+1]=1
		UC2z[nbands+2*i+1][nbands+2*i]=1





	blockzeros=np.full((nbands,nbands),0)
	return UC2z
# C3 symmetry operator in one spin flavor
def C3full():
	nbands=np.size(qxvalstot)*2
	UC3=np.full((nbands,nbands),0.+0.j)
	UC3minus=np.full((nbands,nbands),0.+0.j)
	for i in range(np.size(qxvalstot)):
		kxprime,kyprime=C3rotate(qxvalstot[i],qyvalstot[i])


		x=np.where(np.round(qxvalstot/ktheta,3)==round(kxprime/ktheta, 3))
		y=np.where(np.round(qyvalstot/ktheta,3)==round(kyprime/ktheta, 3))


		for j in range( np.size(np.intersect1d(x, y))):
			if qvalslayer[i]==qvalslayer[np.intersect1d(x, y)[j]]:
				c3index=np.intersect1d(x, y)[j]

		UC3[2*c3index][2*i]=np.exp(-1j*2*np.pi/3)
		UC3[2*c3index+1][2*i+1]=np.exp(1j*2*np.pi/3)

	for i in range(np.size(qxvalstotminus)):
		kxprime,kyprime=C3rotate(qxvalstotminus[i],qyvalstotminus[i])


		x=np.where(np.round(qxvalstotminus/ktheta,3)==round(kxprime/ktheta, 3))
		y=np.where(np.round(qyvalstotminus/ktheta,3)==round(kyprime/ktheta, 3))


		for j in range( np.size(np.intersect1d(x, y))):
			if qvalslayerminus[i]==qvalslayerminus[np.intersect1d(x, y)[j]]:
				c3index=np.intersect1d(x, y)[j]

		UC3minus[2*c3index][2*i]=np.exp(1j*2*np.pi/3)
		UC3minus[2*c3index+1][2*i+1]=np.exp(-1j*2*np.pi/3)

	blockzeros=np.full((nbands,nbands),0)
	UC3Full=np.block([[UC3,blockzeros],[blockzeros,UC3minus]])



	return UC3Full
# C3 symmetry operator with Spin
def C3Spin():
	nbands=np.size(qxvalstot)*2
	UC3=np.full((nbands,nbands),0.+0.j)
	UC3minus=np.full((nbands,nbands),0.+0.j)
	for i in range(np.size(qxvalstot)):
		kxprime,kyprime=C3rotate(qxvalstot[i],qyvalstot[i])


		x=np.where(np.round(qxvalstot/ktheta,3)==round(kxprime/ktheta, 3))
		y=np.where(np.round(qyvalstot/ktheta,3)==round(kyprime/ktheta, 3))


		for j in range( np.size(np.intersect1d(x, y))):
			if qvalslayer[i]==qvalslayer[np.intersect1d(x, y)[j]]:
				c3index=np.intersect1d(x, y)[j]

		UC3[2*c3index][2*i]=np.exp(-1j*2*np.pi/3)
		UC3[2*c3index+1][2*i+1]=np.exp(1j*2*np.pi/3)

	for i in range(np.size(qxvalstotminus)):
		kxprime,kyprime=C3rotate(qxvalstotminus[i],qyvalstotminus[i])


		x=np.where(np.round(qxvalstotminus/ktheta,3)==round(kxprime/ktheta, 3))
		y=np.where(np.round(qyvalstotminus/ktheta,3)==round(kyprime/ktheta, 3))


		for j in range( np.size(np.intersect1d(x, y))):
			if qvalslayerminus[i]==qvalslayerminus[np.intersect1d(x, y)[j]]:
				c3index=np.intersect1d(x, y)[j]

		UC3minus[2*c3index][2*i]=np.exp(1j*2*np.pi/3)
		UC3minus[2*c3index+1][2*i+1]=np.exp(-1j*2*np.pi/3)

	blockzeros=np.full((nbands,nbands),0)
	UC3Full=np.block([[UC3,blockzeros,blockzeros,blockzeros],[blockzeros,UC3minus,blockzeros,blockzeros],[blockzeros,blockzeros,UC3,blockzeros],[blockzeros,blockzeros,blockzeros,UC3minus]])



	return UC3Full

#get ordered eigenvectors for not necessarily hermitian operator (might not be orthonormal if degeneracies)
def eigsystem(H):
    w, v = np.linalg.eig(H)
    idx = np.argsort(w)
    w = w[idx]
    v = v[:,idx]
    return w,v
def eigsystem0(H):
	M = Matrix(H)
	P, D = M.diagonalize()


	w=np.array(D).diagonal()
	wre=[float(re(item1)) for item1 in w]
	wim=[float(im(item1)) for item1 in w]
	w=np.array(wre)+1j*np.array(wim)
	v=np.array(P)
	vre=[[float(re(item2)) for item2 in item1] for item1 in v]
	vim=[[float(im(item2)) for item2 in item1] for item1 in v]
	v=np.array(vre)+1j*np.array(vim)
	idx = np.argsort(w)
	w = w[idx]
	v = v[:,idx]
	return w,v
#get ordered eigenvectors for hermitian operator 
def eighsystem(H):
    w, v = np.linalg.eigh(H)
    idx = np.argsort(w)
    w = w[idx]
    v = v[:,idx]
    return w,v
# Define k independent Moire potential part of Hamiltonian
def MoireMatrix():
	nbands=np.size(qxvalstot)*2
	mat=np.full((nbands,nbands),0.+0.j)
	count=0
	for i in range(np.size(qxvalstot)):
		for j in range(np.size(qxvalstot)):

			if np.abs(ktheta-np.abs(np.sqrt((qxvalstot[i]-qxvalstot[j])**2+(qyvalstot[i]-qyvalstot[j])**2)))<10**(-4)*ktheta:
				if (np.abs(qxvalstot[i]-qxvalstot[j]-q1x)<10**(-3)*ktheta and np.abs(qyvalstot[i]-qyvalstot[j]-q1y)<10**(-3)*ktheta) and (qvalslayer[i]==1 or qvalslayer[j]==1):


					mat[2*i,2*j]=T12_1[0][0]
					mat[2*i,2*j+1]=T12_1[0][1]
					mat[2*i+1,2*j]=T12_1[1][0]
					mat[2*i+1,2*j+1]=T12_1[1][1]

				if  (np.abs(qxvalstot[i]-qxvalstot[j]+q1x)<10**(-3)*ktheta and np.abs(qyvalstot[i]-qyvalstot[j]+q1y)<10**(-3)*ktheta) and (qvalslayer[i]==1 or qvalslayer[j]==1):


					mat[2*i,2*j]=T12_1c[0][0]
					mat[2*i,2*j+1]=T12_1c[0][1]
					mat[2*i+1,2*j]=T12_1c[1][0]
					mat[2*i+1,2*j+1]=T12_1c[1][1]

				if (np.abs(qxvalstot[i]-qxvalstot[j]-q2x)<10**(-3)*ktheta and np.abs(qyvalstot[i]-qyvalstot[j]-q2y)<10**(-3)*ktheta) and (qvalslayer[i]==1 or qvalslayer[j]==1):

					mat[2*i,2*j]=T12_2[0][0]
					mat[2*i,2*j+1]=T12_2[0][1]
					mat[2*i+1,2*j]=T12_2[1][0]
					mat[2*i+1,2*j+1]=T12_2[1][1]

				if (np.abs(qxvalstot[i]-qxvalstot[j]+q2x)<10**(-3)*ktheta and np.abs(qyvalstot[i]-qyvalstot[j]+q2y)<10**(-3)*ktheta) and (qvalslayer[i]==1 or qvalslayer[j]==1):

					mat[2*i,2*j]=T12_2c[0][0]
					mat[2*i,2*j+1]=T12_2c[0][1]
					mat[2*i+1,2*j]=T12_2c[1][0]
					mat[2*i+1,2*j+1]=T12_2c[1][1]

				if (np.abs(qxvalstot[i]-qxvalstot[j]-q3x)<10**(-3)*ktheta and np.abs(qyvalstot[i]-qyvalstot[j]-q3y)<10**(-3)*ktheta)  and (qvalslayer[i]==1 or qvalslayer[j]==1):

					mat[2*i,2*j]=T12_3[0][0]
					mat[2*i,2*j+1]=T12_3[0][1]
					mat[2*i+1,2*j]=T12_3[1][0]
					mat[2*i+1,2*j+1]=T12_3[1][1]

				if  (np.abs(qxvalstot[i]-qxvalstot[j]+q3x)<10**(-3)*ktheta and np.abs(qyvalstot[i]-qyvalstot[j]+q3y)<10**(-3)*ktheta) and (qvalslayer[i]==1 or qvalslayer[j]==1):

					mat[2*i,2*j]=T12_3c[0][0]
					mat[2*i,2*j+1]=T12_3c[0][1]
					mat[2*i+1,2*j]=T12_3c[1][0]
					mat[2*i+1,2*j+1]=T12_3c[1][1]

				if (np.abs(qxvalstot[i]-qxvalstot[j]-q1x)<10**(-3)*ktheta and np.abs(qyvalstot[i]-qyvalstot[j]-q1y)<10**(-3)*ktheta) and (qvalslayer[i]==3 or qvalslayer[j]==3):

					mat[2*i,2*j]=T23_1c[0][0]
					mat[2*i,2*j+1]=T23_1c[0][1]
					mat[2*i+1,2*j]=T23_1c[1][0]
					mat[2*i+1,2*j+1]=T23_1c[1][1]

				if (np.abs(qxvalstot[i]-qxvalstot[j]+q1x)<10**(-3)*ktheta and np.abs(qyvalstot[i]-qyvalstot[j]+q1y)<10**(-3)*ktheta) and (qvalslayer[i]==3 or qvalslayer[j]==3):

					mat[2*i,2*j]=T23_1[0][0]
					mat[2*i,2*j+1]=T23_1[0][1]
					mat[2*i+1,2*j]=T23_1[1][0]
					mat[2*i+1,2*j+1]=T23_1[1][1]

				if (np.abs(qxvalstot[i]-qxvalstot[j]-q2x)<10**(-3)*ktheta and np.abs(qyvalstot[i]-qyvalstot[j]-q2y)<10**(-3)*ktheta) and (qvalslayer[i]==3 or qvalslayer[j]==3):


					mat[2*i,2*j]=T23_2c[0][0]
					mat[2*i,2*j+1]=T23_2c[0][1]
					mat[2*i+1,2*j]=T23_2c[1][0]
					mat[2*i+1,2*j+1]=T23_2c[1][1]

				if (np.abs(qxvalstot[i]-qxvalstot[j]+q2x)<10**(-3)*ktheta and np.abs(qyvalstot[i]-qyvalstot[j]+q2y)<10**(-3)*ktheta) and (qvalslayer[i]==3 or qvalslayer[j]==3):


					mat[2*i,2*j]=T23_2[0][0]
					mat[2*i,2*j+1]=T23_2[0][1]
					mat[2*i+1,2*j]=T23_2[1][0]
					mat[2*i+1,2*j+1]=T23_2[1][1]

				if (np.abs(qxvalstot[i]-qxvalstot[j]-q3x)<10**(-3)*ktheta and np.abs(qyvalstot[i]-qyvalstot[j]-q3y)<10**(-3)*ktheta) and (qvalslayer[i]==3 or qvalslayer[j]==3):

					mat[2*i,2*j]=T23_3c[0][0]
					mat[2*i,2*j+1]=T23_3c[0][1]
					mat[2*i+1,2*j]=T23_3c[1][0]
					mat[2*i+1,2*j+1]=T23_3c[1][1]

				if (np.abs(qxvalstot[i]-qxvalstot[j]+q3x)<10**(-3)*ktheta and np.abs(qyvalstot[i]-qyvalstot[j]+q3y)<10**(-3)*ktheta) and (qvalslayer[i]==3 or qvalslayer[j]==3):

					mat[2*i,2*j]=T23_3[0][0]
					mat[2*i,2*j+1]=T23_3[0][1]
					mat[2*i+1,2*j]=T23_3[1][0]
					mat[2*i+1,2*j+1]=T23_3[1][1]
	return mat
# 2x2 matrix for just 1 dirac cone
def KineticTerm(kx,ky,theta):
	sigma1=np.cos(theta/4)*Id2+1j*pauliz*np.sin(theta/4)
	sigma1c=np.cos(theta/4)*Id2-1j*pauliz*np.sin(theta/4)
	sigmatheta=vF*sigma1.dot(-pauliy*ky+paulix*kx).dot(sigma1c)
	return sigmatheta
# Define full TTG matrix
def Hplus(kx,ky):
	theta=0
	nbands=np.size(qxvalstot)*2
	kinetic=np.full((nbands,nbands),0.+0.j)

	
	for i in range(np.size(qxvalstot)):


		if qvalslayer[i]==1:	
			kinetic[2*i,2*i]=KineticTerm(kx-qxvalstot[i],ky-qyvalstot[i],theta)[0][0]+D
			kinetic[2*i,2*i+1]=KineticTerm(kx-qxvalstot[i],ky-qyvalstot[i],theta)[0][1]
			kinetic[2*i+1,2*i]=KineticTerm(kx-qxvalstot[i],ky-qyvalstot[i],theta)[1][0]
			kinetic[2*i+1,2*i+1]=KineticTerm(kx-qxvalstot[i],ky-qyvalstot[i],theta)[1][1]+D

		elif qvalslayer[i]==2:	
			kinetic[2*i,2*i]=KineticTerm(kx-qxvalstot[i],ky-qyvalstot[i],-theta)[0][0]
			kinetic[2*i,2*i+1]=KineticTerm(kx-qxvalstot[i],ky-qyvalstot[i],-theta)[0][1]
			kinetic[2*i+1,2*i]=KineticTerm(kx-qxvalstot[i],ky-qyvalstot[i],-theta)[1][0]
			kinetic[2*i+1,2*i+1]=KineticTerm(kx-qxvalstot[i],ky-qyvalstot[i],-theta)[1][1]
		elif qvalslayer[i]==3:	
			kinetic[2*i,2*i]=KineticTerm(kx-qxvalstot[i],ky-qyvalstot[i],theta)[0][0]-D
			kinetic[2*i,2*i+1]=KineticTerm(kx-qxvalstot[i],ky-qyvalstot[i],theta)[0][1]
			kinetic[2*i+1,2*i]=KineticTerm(kx-qxvalstot[i],ky-qyvalstot[i],theta)[1][0]
			kinetic[2*i+1,2*i+1]=KineticTerm(kx-qxvalstot[i],ky-qyvalstot[i],theta)[1][1]-D
		else:
			kinetic[2*i,2*i]=KineticTerm(kx-qxvalstot[i],ky-qyvalstot[i],theta)[0][0]
			kinetic[2*i,2*i+1]=KineticTerm(kx-qxvalstot[i],ky-qyvalstot[i],theta)[0][1]
			kinetic[2*i+1,2*i]=KineticTerm(kx-qxvalstot[i],ky-qyvalstot[i],theta)[1][0]
			kinetic[2*i+1,2*i+1]=KineticTerm(kx-qxvalstot[i],ky-qyvalstot[i],theta)[1][1]

	Hfull=kinetic+moirepotential

	return Hfull
# Define time reversed TTG copy for other valley
def Hminus(kx,ky):
	nbands=np.size(qxvalstot)*2
	blockzeros=np.full((nbands,nbands),0)
	return T0().dot(np.block([[np.conjugate(Hplus(-kx,-ky)),blockzeros],[blockzeros,blockzeros]])).dot(T0())[nbands:2*nbands,nbands:2*nbands]
# Define full H per spin flavor
def H(kx,ky):
	nbands=np.size(qxvalstot)*2
	blockzeros=np.full((nbands,nbands),0)
	return np.block([[Hplus(kx,ky),blockzeros],[blockzeros,blockzeros]])+T0().dot(np.block([[np.conjugate(Hplus(-kx,-ky)),blockzeros],[blockzeros,blockzeros]])).dot(T0())
# Mz reflection operator
def MZ():
	nbands=np.size(qxvalstot)*2
	U=np.full((nbands,nbands),0.+0.j)
	Uminus=np.full((nbands,nbands),0.+0.j)
	for i in range(np.size(qxvalstot)):
		for j in range(np.size(qxvalstot)):

			if qvalslayer[i]==2 and qvalslayer[j]==2 and qxvalstot[i]==qxvalstot[j] and qyvalstot[i]==qyvalstot[j]:
				U[2*i][2*j]=1
				U[2*i+1][2*j+1]=1


			if qvalslayer[i]==1 and qvalslayer[j]==3 and qxvalstot[i]==qxvalstot[j] and qyvalstot[i]==qyvalstot[j]:
				U[2*i][2*j]=1
				U[2*i+1][2*j+1]=1
			if qvalslayer[i]==3 and qvalslayer[j]==1 and qxvalstot[i]==qxvalstot[j] and qyvalstot[i]==qyvalstot[j]:
				U[2*i][2*j]=1
				U[2*i+1][2*j+1]=1

	for i in range(np.size(qxvalstotminus)):
		for j in range(np.size(qxvalstotminus)):


			if qvalslayerminus[i]==2 and qvalslayerminus[j]==2 and qxvalstotminus[i]==qxvalstotminus[j] and qyvalstotminus[i]==qyvalstotminus[j]:
				Uminus[2*i][2*j]=1
				Uminus[2*i+1][2*j+1]=1


			if qvalslayerminus[i]==1 and qvalslayerminus[j]==3 and qxvalstotminus[i]==qxvalstotminus[j] and qyvalstotminus[i]==qyvalstotminus[j]:
				Uminus[2*i][2*j]=1
				Uminus[2*i+1][2*j+1]=1
			if qvalslayerminus[i]==3 and qvalslayerminus[j]==1 and qxvalstotminus[i]==qxvalstotminus[j] and qyvalstotminus[i]==qyvalstotminus[j]:
				Uminus[2*i][2*j]=1
				Uminus[2*i+1][2*j+1]=1


	blockzeros=np.full((nbands,nbands),0)
	UFull=np.block([[Uminus,blockzeros,blockzeros,blockzeros],[blockzeros,Uminus,blockzeros,blockzeros],[blockzeros,blockzeros,Uminus,blockzeros],[blockzeros,blockzeros,blockzeros,Uminus]])

	return UFull
# PH*C2z symmetry operator
def PHC2z():
	nbands=np.size(qxvalstot)*2
	Uph=np.full((nbands*2,nbands*2),0.+0.j)
	for i in range(np.size(qxvalstot)):
		kxprime,kyprime=qxvalstot[i],qyvalstot[i]


		x=np.where(np.round(qxvalstotminus/ktheta,3)==round(kxprime/ktheta, 3))
		y=np.where(np.round(qyvalstotminus/ktheta,3)==round(kyprime/ktheta, 3))

		phindex=0

		intersect=np.intersect1d(x, y)
		# 2 and 1 coupling only

		for j in range( np.size(intersect)):
			if( qvalslayer[i]==2 and qvalslayerminus[np.intersect1d(x, y)[j]]==1) or (qvalslayer[i]==2 and qvalslayerminus[np.intersect1d(x, y)[j]]==3) :
				phindex=np.intersect1d(x, y)[j]

				Uph[2*i][nbands+2*phindex+1]=-1j*1j*np.sqrt(2)
				Uph[2*i+1][nbands+2*phindex]=-1j*1j*np.sqrt(2)
				Uph[nbands+2*phindex][2*i+1]=1j*1j*np.sqrt(2)
				Uph[nbands+2*phindex+1][2*i]=1j*1j*np.sqrt(2)
			if (qvalslayer[i]==1 and qvalslayerminus[np.intersect1d(x, y)[j]]==2) or (qvalslayer[i]==3 and qvalslayerminus[np.intersect1d(x, y)[j]]==2):
				phindex=np.intersect1d(x, y)[j]

				Uph[2*i][nbands+2*phindex+1]=1j*1j*np.sqrt(2)
				Uph[2*i+1][nbands+2*phindex]=1j*1j*np.sqrt(2)
				Uph[nbands+2*phindex][2*i+1]=-1j*1j*np.sqrt(2)
				Uph[nbands+2*phindex+1][2*i]=-1j*1j*np.sqrt(2)

		x=np.where(np.round(qxvalstot/ktheta,3)==round(kxprime/ktheta, 3))
		y=np.where(np.round(qyvalstot/ktheta,3)==round(kyprime/ktheta, 3))

		phindex=0

		# 3 and 1 coupling only
		intersect=np.intersect1d(x, y)
		for j in range( np.size(intersect)):
			if qvalslayer[i]==1 and qvalslayer[np.intersect1d(x, y)[j]]==3:
				phindex=np.intersect1d(x, y)[j]
				Uph[2*i][2*phindex]=-1j
				Uph[2*i+1][2*phindex+1]=1j

			if qvalslayer[i]==3 and qvalslayer[np.intersect1d(x, y)[j]]==1:
				phindex=np.intersect1d(x, y)[j]
				Uph[2*i][2*phindex]=-1j
				Uph[2*i+1][2*phindex+1]=1j


			if qvalslayer[i]==1 and qvalslayer[np.intersect1d(x, y)[j]]==1:
				phindex=np.intersect1d(x, y)[j]
				Uph[2*i][2*phindex]=1j
				Uph[2*i+1][2*phindex+1]=-1j



			if qvalslayer[i]==3 and qvalslayer[np.intersect1d(x, y)[j]]==3:
				phindex=np.intersect1d(x, y)[j]
				Uph[2*i][2*phindex]=1j
				Uph[2*i+1][2*phindex+1]=-1j

		kxprime,kyprime=qxvalstotminus[i],qyvalstotminus[i]

		x=np.where(np.round(qxvalstotminus/ktheta,3)==round(kxprime/ktheta, 3))
		y=np.where(np.round(qyvalstotminus/ktheta,3)==round(kyprime/ktheta, 3))

		phindex=0

		# 3 and 1 coupling only
		intersect=np.intersect1d(x, y)
		for j in range( np.size(intersect)):
			if qvalslayerminus[i]==1 and qvalslayerminus[np.intersect1d(x, y)[j]]==3:
				phindex=np.intersect1d(x, y)[j]
				Uph[nbands+2*i][nbands+2*phindex]=1j
				Uph[nbands+2*i+1][nbands+2*phindex+1]=-1j

			if qvalslayerminus[i]==3 and qvalslayerminus[np.intersect1d(x, y)[j]]==1:
				phindex=np.intersect1d(x, y)[j]
				Uph[nbands+2*i][nbands+2*phindex]=1j
				Uph[nbands+2*i+1][nbands+2*phindex+1]=-1j


			if qvalslayerminus[i]==1 and qvalslayerminus[np.intersect1d(x, y)[j]]==1:
				phindex=np.intersect1d(x, y)[j]
				Uph[nbands+2*i][nbands+2*phindex]=-1j
				Uph[nbands+2*i+1][nbands+2*phindex+1]=1j



			if qvalslayerminus[i]==3 and qvalslayerminus[np.intersect1d(x, y)[j]]==3:
				phindex=np.intersect1d(x, y)[j]
				Uph[nbands+2*i][nbands+2*phindex]=-1j
				Uph[nbands+2*i+1][nbands+2*phindex+1]=1j


	Uph=Uph*1/2




	blockzeros=np.full((nbands,nbands),0)

	return Uph
# Chiral symmetry operator
def Chiral():

	nbands=np.size(qxvalstot)*2
	Uchiral=np.full((nbands*2,nbands*2),0.+0.j)
	for i in range(np.size(qxvalstot)):
		kxprime=qxvalstot[i]
		kyprime=qyvalstot[i]


		x=np.where(np.round(qxvalstot/ktheta,3)==round(kxprime/ktheta, 3))
		y=np.where(np.round(qyvalstot/ktheta,3)==round(kyprime/ktheta, 3))


		intersect=np.intersect1d(x, y)

		for j in range( np.size(intersect)):
			if qvalslayer[i]!=qvalslayer[np.intersect1d(x, y)[j]]:

				index2=np.intersect1d(x, y)[j]

				Uchiral[2*i][2*index2]=1
				Uchiral[2*i+1][2*index2+1]=-1

			if qvalslayer[i]==qvalslayer[np.intersect1d(x, y)[j]] and qvalslayer[i]==2:

				
				index2=np.intersect1d(x, y)[j]

				Uchiral[2*i][2*index2]=1
				Uchiral[2*i+1][2*index2+1]=-1

	for i in range(np.size(qxvalstotminus)):
		kxprime=qxvalstotminus[i]
		kyprime=qyvalstotminus[i]


		x=np.where(np.round(qxvalstotminus/ktheta,3)==round(kxprime/ktheta, 3))
		y=np.where(np.round(qyvalstotminus/ktheta,3)==round(kyprime/ktheta, 3))


		intersect=np.intersect1d(x, y)

		for j in range( np.size(intersect)):
			if qvalslayerminus[i]!=qvalslayerminus[np.intersect1d(x, y)[j]]:
				index2=np.intersect1d(x, y)[j]

				Uchiral[2*i+nbands][2*index2+nbands]=1
				Uchiral[2*i+1+nbands][2*index2+1+nbands]=-1

			if qvalslayerminus[i]==qvalslayerminus[np.intersect1d(x, y)[j]] and qvalslayerminus[i]==2:
				index2=np.intersect1d(x, y)[j]

				Uchiral[2*i+nbands][2*index2+nbands]=1
				Uchiral[2*i+1+nbands][2*index2+1+nbands]=-1




	blockzeros=np.full((2*nbands,2*nbands),0)
	UchiralFull=np.block([[Uchiral,blockzeros],[blockzeros,Uchiral]])
	return Uchiral
# T symmetry operator in one spin flavor
def T0():
	nbands=np.size(qxvalstot)*2
	UT=np.full((nbands*2,nbands*2),0.+0.j)
	for i in range(np.size(qxvalstot)):
		kxprime,kyprime=C2zrotate(qxvalstot[i],qyvalstot[i])


		x=np.where(np.round(qxvalstotminus/ktheta,3)==round(kxprime/ktheta, 3))
		y=np.where(np.round(qyvalstotminus/ktheta,3)==round(kyprime/ktheta, 3))

		c2zindex=0

		intersect=np.intersect1d(x, y)

		for j in range( np.size(intersect)):
			if qvalslayer[i]==qvalslayerminus[np.intersect1d(x, y)[j]]:
				c2zindex=np.intersect1d(x, y)[j]


		UT[2*i][nbands+2*c2zindex]=1
		UT[2*i+1][nbands+2*c2zindex+1]=1
		UT[nbands+2*c2zindex][2*i]=1
		UT[nbands+2*c2zindex+1][2*i+1]=1



	return UT

############################################################################################################################################
moirepotential=MoireMatrix()
# given number of pt in brillouin zone i, return momentum (kx,ky)
kxvals=[]
kyvals=[]
kxvalsnew=[]
kyvalsnew=[]

# Make first half of k grid
for i in range(nptsd):
	for j in range(nptsd+i):
		if i==nptsd-1 and (j==0 or j==2*nptsd-2):
			continue
		if i==nptsd-1 and j==0:
			continue


		kxvals.append(-vy-by*i-c1y*j)
		kyvals.append(vx+bx*i+c1x*j)

# Make second half of k grid
for i in range(nptsd-1):
	for j in range(2*nptsd-2-i):
		if j==0 or j==2*nptsd-3-i:
			continue
		if i == nptsd-2:
			continue

		kxvals.append(-vy-(nptsd-1)*by-c2y*(1+i)-c1y*j)
		kyvals.append(vx+(nptsd-1)*bx+c2x*(1+i)+c1x*j)

# The below (line 1080-1241) is for plotting purposes only, to check the G-grid was made successfully. Ignore for HF calculations
kxc61=[]
kyc61=[]
kxc62=[]
kyc62=[]
kxc63=[]
kyc63=[]
kxc64=[]
kyc64=[]
kxc65=[]
kyc65=[]
kxc66=[]
kyc66=[]

for j in range(np.size(kxvals)):



	if kxvals[j]/ktheta<0-10**(-3) and kxvals[j]/ktheta*q3y/q3x-kyvals[j]/ktheta<10**(-3):
		kxc61.append(kxvals[j])
		kyc61.append(kyvals[j])

	elif kxvals[j]/ktheta>0-10**(-3) and -kxvals[j]/ktheta*q3y/q3x-kyvals[j]/ktheta<-10**(-3):
		kxc62.append(kxvals[j])
		kyc62.append(kyvals[j])

	elif kxvals[j]/ktheta>0+10**(-3) and -kxvals[j]/ktheta*q3y/q3x-kyvals[j]/ktheta>-10**(-3) and kxvals[j]/ktheta*q3y/q3x-kyvals[j]/ktheta<-10**(-3):
		kxc63.append(kxvals[j])
		kyc63.append(kyvals[j])

	elif kxvals[j]/ktheta>0+10**(-3) and kxvals[j]/ktheta*q3y/q3x-kyvals[j]/ktheta>-10**(-3):
		kxc64.append(kxvals[j])
		kyc64.append(kyvals[j])

	elif kxvals[j]/ktheta<0+10**(-3) and -kxvals[j]/ktheta*q3y/q3x-kyvals[j]/ktheta>10**(-3):
		kxc65.append(kxvals[j])
		kyc65.append(kyvals[j])
	elif kxvals[j]/ktheta<0-10**(-3) and -kxvals[j]/ktheta*q3y/q3x-kyvals[j]/ktheta<10**(-3) and kxvals[j]/ktheta*q3y/q3x-kyvals[j]/ktheta>10**(-3):
		kxc66.append(kxvals[j])
		kyc66.append(kyvals[j])


kxshell1=[]
kyshell1=[]
for i in range(np.size(shell1momentax)):
	for j in range(np.size(kxvals)):

		if np.sqrt((shell1momentax[i]+kxvals[j])**2+(shell1momentay[i]+kyvals[j])**2)/ktheta >2-.1 and (np.abs(kxvals[j]+b1x)<10**(-3)*ktheta or np.abs(kyvals[j]-kxvals[j]*q2y/q2x-ktheta)<10**(-3) or np.abs(kyvals[j]+kxvals[j]*q2y/q2x-ktheta)<10**(-3)):
			continue

		if (kxvals[j]+shell1momentax[i])/ktheta<0-10**(-3) and (kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta<10**(-3):
			kxc61.append((kxvals[j]+shell1momentax[i]))
			kyc61.append((kyvals[j]+shell1momentay[i]))

		elif (kxvals[j]+shell1momentax[i])/ktheta>0-10**(-3) and -(kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta<-10**(-3):
			kxc62.append((kxvals[j]+shell1momentax[i]))
			kyc62.append((kyvals[j]+shell1momentay[i]))

		elif (kxvals[j]+shell1momentax[i])/ktheta>0+10**(-3) and -(kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta>-10**(-3) and (kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta<-10**(-3):
			kxc63.append((kxvals[j]+shell1momentax[i]))
			kyc63.append((kyvals[j]+shell1momentay[i]))

		elif (kxvals[j]+shell1momentax[i])/ktheta>0+10**(-3) and (kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta>-10**(-3):
			kxc64.append((kxvals[j]+shell1momentax[i]))
			kyc64.append((kyvals[j]+shell1momentay[i]))

		elif (kxvals[j]+shell1momentax[i])/ktheta<0+10**(-3) and -(kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta>10**(-3):
			kxc65.append((kxvals[j]+shell1momentax[i]))
			kyc65.append((kyvals[j]+shell1momentay[i]))
		elif (kxvals[j]+shell1momentax[i])/ktheta<0-10**(-3) and -(kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta<10**(-3) and (kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta>10**(-3):
			kxc66.append((kxvals[j]+shell1momentax[i]))
			kyc66.append((kyvals[j]+shell1momentay[i]))
for i in range(np.size(shell2momentax)):
	for j in range(np.size(kxvals)):

		

		if (kxvals[j]+shell2momentax[i])/ktheta<0-10**(-3) and (kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta<10**(-3):
			kxc61.append((kxvals[j]+shell2momentax[i]))
			kyc61.append((kyvals[j]+shell2momentay[i]))

		elif (kxvals[j]+shell2momentax[i])/ktheta>0-10**(-3) and -(kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta<-10**(-3):
			kxc62.append((kxvals[j]+shell2momentax[i]))
			kyc62.append((kyvals[j]+shell2momentay[i]))

		elif (kxvals[j]+shell2momentax[i])/ktheta>0+10**(-3) and -(kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta>-10**(-3) and (kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta<-10**(-3):
			kxc63.append((kxvals[j]+shell2momentax[i]))
			kyc63.append((kyvals[j]+shell2momentay[i]))

		elif (kxvals[j]+shell2momentax[i])/ktheta>0+10**(-3) and (kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta>-10**(-3):
			kxc64.append((kxvals[j]+shell2momentax[i]))
			kyc64.append((kyvals[j]+shell2momentay[i]))

		elif (kxvals[j]+shell2momentax[i])/ktheta<0+10**(-3) and -(kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta>10**(-3):
			kxc65.append((kxvals[j]+shell2momentax[i]))
			kyc65.append((kyvals[j]+shell2momentay[i]))
		elif (kxvals[j]+shell2momentax[i])/ktheta<0-10**(-3) and -(kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta<10**(-3) and (kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta>10**(-3):
			kxc66.append((kxvals[j]+shell2momentax[i]))
			kyc66.append((kyvals[j]+shell2momentay[i]))

kqtestx=[]
kqtesty=[]
for i in range(npts):
	for j in range(npts):
		for k in range(np.size(shell1momentax)):
			kqtestx.append(kxvals[i]+kxvals[j])
			kqtesty.append(kyvals[i]+kyvals[j])
			if (np.abs(kxvals[i]+b1x)<10**(-3)*ktheta or np.abs(kyvals[i]-kxvals[i]*q2y/q2x-ktheta)<10**(-3) or np.abs(kyvals[i]+kxvals[i]*q2y/q2x-ktheta)<10**(-3)):
				continue
			if np.sqrt((shell1momentax[k]+kxvals[j])**2+(shell1momentay[k]+kyvals[j])**2)/ktheta >2-.1 and (np.abs(kxvals[j]+b1x)<10**(-3)*ktheta or np.abs(kyvals[j]-kxvals[j]*q2y/q2x-ktheta)<10**(-3) or np.abs(kyvals[j]+kxvals[j]*q2y/q2x-ktheta)<10**(-3)):
					continue
			else:
				kqtestx.append(kxvals[i]+kxvals[j]+shell1momentax[k])
				kqtesty.append(kyvals[i]+kyvals[j]+shell1momentay[k])


kqbvector=[]
kqindex=[]
kxsvalsnew=[]
kysvalsnew=[]


# construct scattering momentum matrix given two pts in reduced BZ
for i in range(npts):
	kqindexx=[]
	kqbvectorx=[]
	for j in range(npts):
		#get k and q from i and j
		kx=kxvals[i]
		ky=kyvals[i]
		qx=kxvals[j]
		qy=kyvals[j]
		scatx=(kx+qx)
		scaty=(ky+qy)


		bvectorx=0
		bvectory=0

		if(-scaty/scatx>=-1/np.sqrt(3)+10**(-3) and -scaty/scatx<1/np.sqrt(3)-10**(-3) and -scatx>np.sqrt(3)/2*ktheta+10**(-3)):
			scatx=scatx+b1x+b2x

			scatxprime=scatx
			scatyprime=scaty
			kxsvalsnew.append(scatxprime)
			kysvalsnew.append(scatyprime)

			bvectorx=-b1x-b2x
			bvectory=-b1y-b2y
		if(-scaty/scatx>=-1/np.sqrt(3)+10**(-3) and -scaty/scatx<1/np.sqrt(3)+10**(-3) and -scatx>np.sqrt(3)/2*ktheta+10**(-3)):
			scatx=scatx+b1x+b2x

			scatxprime=scatx
			scatyprime=scaty
			kxsvalsnew.append(scatxprime)
			kysvalsnew.append(scatyprime)

			bvectorx=-b1x-b2x
			bvectory=-b1y-b2y


		
		if(-scaty/scatx>-1/np.sqrt(3)+10**(-3) and -scaty/scatx<1/np.sqrt(3)-10**(-3) and -scatx<-np.sqrt(3)/2*ktheta+10**(-3) ):
			scatx=scatx-b1x-b2x

			scatxprime=scatx
			scatyprime=scaty

			bvectorx=b1x+b2x
			bvectory=b1y+b2y
			kxsvalsnew.append(scatxprime)
			kysvalsnew.append(scatyprime)
		if(-scaty/scatx>-1/np.sqrt(3)-10**(-3) and -scaty/scatx<1/np.sqrt(3)-10**(-3) and -scatx<-np.sqrt(3)/2*ktheta+10**(-3) ):
			scatx=scatx-b1x-b2x

			scatxprime=scatx
			scatyprime=scaty

			bvectorx=b1x+b2x
			bvectory=b1y+b2y
			kxsvalsnew.append(scatxprime)
			kysvalsnew.append(scatyprime)




		if(-scatx/scaty>=-np.sqrt(3)+10**(-3) and -scatx<-10**(-3) and scaty > ktheta-scatx/np.sqrt(3)+10**(-3)):
			scatx=scatx-b1x
			scaty=scaty-b1y

			bvectorx=b1x
			bvectory=b1y

			scatxprime=scatx
			scatyprime=scaty

		if(-scatx/scaty>=-np.sqrt(3)+10**(-3) and -scatx<+10**(-3) and scaty > ktheta-scatx/np.sqrt(3)+10**(-3)):
			scatx=scatx-b1x
			scaty=scaty-b1y

			bvectorx=b1x
			bvectory=b1y

			scatxprime=scatx
			scatyprime=scaty



		if(-scatx/scaty>=-np.sqrt(3)+10**(-3) and -scatx>+10**(-3) and scaty <= -ktheta-scatx/np.sqrt(3)+10**(-3)):
			scatx=scatx+b1x
			scaty=scaty+b1y

			bvectorx=-b1x
			bvectory=-b1y
			# scatx=10**(-8)
			# scaty=10**(-8)
			scatxprime=scatx
			scatyprime=scaty

		if(-scatx/scaty>=-np.sqrt(3)-10**(-3) and -scatx>+10**(-3) and scaty <= -ktheta-scatx/np.sqrt(3)+10**(-3)):
			scatx=scatx+b1x
			scaty=scaty+b1y


			bvectorx=-b1x
			bvectory=-b1y
			scatxprime=scatx
			scatyprime=scaty

		

		if(-scatx/scaty<np.sqrt(3)-10**(-3) and -scatx>=10**(-3) and scaty >= ktheta+scatx/np.sqrt(3)+10**(-3)):
			scatx=scatx+b2x
			scaty=scaty+b2y

			bvectorx=-b2x
			bvectory=-b2y

			scatxprime=scatx
			scatyprime=scaty

		if(-scatx/scaty<np.sqrt(3)-10**(-3) and -scatx>=10**(-3) and scaty >= ktheta+scatx/np.sqrt(3)+10**(-3)):
			scatx=scatx+b2x
			scaty=scaty+b2y

			bvectorx=-b2x
			bvectory=-b2y

			scatxprime=scatx
			scatyprime=scaty


		if(-scatx/scaty<=np.sqrt(3)-10**(-3) and -scatx<-10**(-3) and scaty <= -ktheta+scatx/np.sqrt(3)+10**(-3)):
			scatx=scatx-b2x
			scaty=scaty-b2y

			bvectorx=b2x
			bvectory=b2y

			scatxprime=scatx
			scatyprime=scaty

		if(-scatx/scaty<=np.sqrt(3)+10**(-3) and -scatx<+10**(-3) and scaty <= -ktheta+scatx/np.sqrt(3)+10**(-3)):
			scatx=scatx-b2x
			scaty=scaty-b2y

			bvectorx=b2x
			bvectory=b2y

			scatxprime=scatx
			scatyprime=scaty




		x=np.where(np.round(kxvals,5)==round(scatx, 5))
		y=np.where(np.round(kyvals,5)==round(scaty, 5))


		xvecs=np.where(np.round(shell1momentax/ktheta,5)==round(bvectorx/ktheta, 5))
		yvecs=np.where(np.round(shell1momentay/ktheta,5)==round(bvectory/ktheta, 5))

		intersectvecs=np.intersect1d(xvecs, yvecs)

		if np.size(intersectvecs)==0:
			kqbvectorx.append(-1)
		else:
			kqbvectorx.append(intersectvecs[0])


		kxvalsnew.append(scatx)
		kyvalsnew.append(scaty)

		index=np.intersect1d(x, y)[0]


		kqindexx.append(np.intersect1d(x, y)[0])
	kqindex.append(kqindexx)
	kqbvector.append(kqbvectorx)
np.save('kqindexrepo',kqindex)
np.save('kqbvectorrepo',kqbvector)

kxsvalsnew=[]
kysvalsnew=[]
kc3index=[]

# construct mapping between two momenta related by a C3 rotation
for i in range(npts):
	kx=kxvals[i]
	ky=kyvals[i]

	kxprime,kyprime=C3rotate(kx,ky)


	if(-kyprime/kxprime>=-1/np.sqrt(3)+10**(-3) and -kyprime/kxprime<1/np.sqrt(3)-10**(-3) and -kxprime>np.sqrt(3)/2*ktheta+10**(-3)):
		kxprime=kxprime+b1x+b2x

		kxprimeprime=kxprime
		kyprimeprime=kyprime
		kxsvalsnew.append(kxprimeprime)
		kysvalsnew.append(kyprimeprime)
	if(-kyprime/kxprime>=-1/np.sqrt(3)+10**(-3) and -kyprime/kxprime<1/np.sqrt(3)+10**(-3) and -kxprime>np.sqrt(3)/2*ktheta+10**(-3)):
		kxprime=kxprime+b1x+b2x

		kxprimeprime=kxprime
		kyprimeprime=kyprime
		kxsvalsnew.append(kxprimeprime)
		kysvalsnew.append(kyprimeprime)


	
	if(-kyprime/kxprime>-1/np.sqrt(3)+10**(-3) and -kyprime/kxprime<1/np.sqrt(3)-10**(-3) and -kxprime<-np.sqrt(3)/2*ktheta+10**(-3) ):
		kxprime=kxprime-b1x-b2x

		kxprimeprime=kxprime
		kyprimeprime=kyprime
		kxsvalsnew.append(kxprimeprime)
		kysvalsnew.append(kyprimeprime)
	if(-kyprime/kxprime>-1/np.sqrt(3)-10**(-3) and -kyprime/kxprime<1/np.sqrt(3)-10**(-3) and -kxprime<-np.sqrt(3)/2*ktheta+10**(-3) ):
		kxprime=kxprime-b1x-b2x

		kxprimeprime=kxprime
		kyprimeprime=kyprime
		kxsvalsnew.append(kxprimeprime)
		kysvalsnew.append(kyprimeprime)




	if(-kxprime/kyprime>=-np.sqrt(3)+10**(-3) and -kxprime<-10**(-3) and kyprime > ktheta-kxprime/np.sqrt(3)+10**(-3)):
		kxprime=kxprime-b1x
		kyprime=kyprime-b1y

		kxprimeprime=kxprime
		kyprimeprime=kyprime

	if(-kxprime/kyprime>=-np.sqrt(3)+10**(-3) and -kxprime<+10**(-3) and kyprime > ktheta-kxprime/np.sqrt(3)+10**(-3)):
		kxprime=kxprime-b1x
		kyprime=kyprime-b1y

		kxprimeprime=kxprime
		kyprimeprime=kyprime


	if(-kxprime/kyprime>=-np.sqrt(3)+10**(-3) and -kxprime>+10**(-3) and kyprime <= -ktheta-kxprime/np.sqrt(3)+10**(-3)):
		kxprime=kxprime+b1x
		kyprime=kyprime+b1y

		kxprimeprime=kxprime
		kyprimeprime=kyprime

	if(-kxprime/kyprime>=-np.sqrt(3)-10**(-3) and -kxprime>+10**(-3) and kyprime <= -ktheta-kxprime/np.sqrt(3)+10**(-3)):
		kxprime=kxprime+b1x
		kyprime=kyprime+b1y

		kxprimeprime=kxprime
		kyprimeprime=kyprime

	

	if(-kxprime/kyprime<np.sqrt(3)-10**(-3) and -kxprime>=10**(-3) and kyprime >= ktheta+kxprime/np.sqrt(3)+10**(-3)):
		kxprime=kxprime+b2x
		kyprime=kyprime+b2y

		kxprimeprime=kxprime
		kyprimeprime=kyprime

	if(-kxprime/kyprime<np.sqrt(3)-10**(-3) and -kxprime>=10**(-3) and kyprime >= ktheta+kxprime/np.sqrt(3)+10**(-3)):
		kxprime=kxprime+b2x
		kyprime=kyprime+b2y

		kxprimeprime=kxprime
		kyprimeprime=kyprime


	if(-kxprime/kyprime<=np.sqrt(3)-10**(-3) and -kxprime<-10**(-3) and kyprime <= -ktheta+kxprime/np.sqrt(3)+10**(-3)):
		kxprime=kxprime-b2x
		kyprime=kyprime-b2y

		kxprimeprime=kxprime
		kyprimeprime=kyprime

	if(-kxprime/kyprime<=np.sqrt(3)+10**(-3) and -kxprime<+10**(-3) and kyprime <= -ktheta+kxprime/np.sqrt(3)+10**(-3)):
		kxprime=kxprime-b2x
		kyprime=kyprime-b2y

		kxprimeprime=kxprime
		kyprimeprime=kyprime


	x=np.where(np.round(kxvals,5)==round(kxprime, 5))
	y=np.where(np.round(kyvals,5)==round(kyprime, 5))



	kxvalsnew.append(kxprime)
	kyvalsnew.append(kyprime)


	kc3index.append(np.intersect1d(x, y)[0])

np.save('kc3indexrepo',kc3index)



kxvalsnew=[]
kyvalsnew=[]
kminusindex=[]

kc2yindex=[]

# construct mapping between two momenta related by a C2z rotation
for i in range(npts):
	#get k and q from i and j
	kx=kxvals[i]
	ky=kyvals[i]
	scatx=-kx
	scaty=-ky

	if(-scaty/scatx>=-1/np.sqrt(3)+10**(-3) and -scaty/scatx<1/np.sqrt(3)-10**(-3) and -scatx>np.sqrt(3)/2*ktheta+10**(-3)):
		scatx=scatx+b1x+b2x

		scatxprime=scatx
		scatyprime=scaty
		kxsvalsnew.append(scatxprime)
		kysvalsnew.append(scatyprime)
	if(-scaty/scatx>=-1/np.sqrt(3)+10**(-3) and -scaty/scatx<1/np.sqrt(3)+10**(-3) and -scatx>np.sqrt(3)/2*ktheta+10**(-3)):
		scatx=scatx+b1x+b2x

		scatxprime=scatx
		scatyprime=scaty
		kxsvalsnew.append(scatxprime)
		kysvalsnew.append(scatyprime)


	
	if(-scaty/scatx>-1/np.sqrt(3)+10**(-3) and -scaty/scatx<1/np.sqrt(3)-10**(-3) and -scatx<-np.sqrt(3)/2*ktheta+10**(-3) ):
		scatx=scatx-b1x-b2x

		scatxprime=scatx
		scatyprime=scaty
		kxsvalsnew.append(scatxprime)
		kysvalsnew.append(scatyprime)
	if(-scaty/scatx>-1/np.sqrt(3)-10**(-3) and -scaty/scatx<1/np.sqrt(3)-10**(-3) and -scatx<-np.sqrt(3)/2*ktheta+10**(-3) ):
		scatx=scatx-b1x-b2x

		scatxprime=scatx
		scatyprime=scaty
		kxsvalsnew.append(scatxprime)
		kysvalsnew.append(scatyprime)




	if(-scatx/scaty>=-np.sqrt(3)+10**(-3) and -scatx<-10**(-3) and scaty > ktheta-scatx/np.sqrt(3)+10**(-3)):
		scatx=scatx-b1x
		scaty=scaty-b1y

		scatxprime=scatx
		scatyprime=scaty

	if(-scatx/scaty>=-np.sqrt(3)+10**(-3) and -scatx<+10**(-3) and scaty > ktheta-scatx/np.sqrt(3)+10**(-3)):
		scatx=scatx-b1x
		scaty=scaty-b1y

		scatxprime=scatx
		scatyprime=scaty



	if(-scatx/scaty>=-np.sqrt(3)+10**(-3) and -scatx>+10**(-3) and scaty <= -ktheta-scatx/np.sqrt(3)+10**(-3)):
		scatx=scatx+b1x
		scaty=scaty+b1y

		scatxprime=scatx
		scatyprime=scaty

	if(-scatx/scaty>=-np.sqrt(3)-10**(-3) and -scatx>+10**(-3) and scaty <= -ktheta-scatx/np.sqrt(3)+10**(-3)):
		scatx=scatx+b1x
		scaty=scaty+b1y

		scatxprime=scatx
		scatyprime=scaty


	

	if(-scatx/scaty<np.sqrt(3)-10**(-3) and -scatx>=10**(-3) and scaty >= ktheta+scatx/np.sqrt(3)+10**(-3)):
		scatx=scatx+b2x
		scaty=scaty+b2y

		scatxprime=scatx
		scatyprime=scaty

	if(-scatx/scaty<np.sqrt(3)-10**(-3) and -scatx>=10**(-3) and scaty >= ktheta+scatx/np.sqrt(3)+10**(-3)):
		scatx=scatx+b2x
		scaty=scaty+b2y

		scatxprime=scatx
		scatyprime=scaty


	if(-scatx/scaty<=np.sqrt(3)-10**(-3) and -scatx<-10**(-3) and scaty <= -ktheta+scatx/np.sqrt(3)+10**(-3)):
		scatx=scatx-b2x
		scaty=scaty-b2y

		scatxprime=scatx
		scatyprime=scaty

	if(-scatx/scaty<=np.sqrt(3)+10**(-3) and -scatx<+10**(-3) and scaty <= -ktheta+scatx/np.sqrt(3)+10**(-3)):
		scatx=scatx-b2x
		scaty=scaty-b2y

		scatxprime=scatx
		scatyprime=scaty



	x=np.where(np.round(kxvals,5)==round(scatx, 5))
	y=np.where(np.round(kyvals,5)==round(scaty, 5))


	kxvalsnew.append(scatx)
	kyvalsnew.append(scaty)


	kminusindex.append(np.intersect1d(x, y)[0])
# REPLACE REFERS TO MY DIRECTORIES, REPLACE WITH OWN PATHS IF USING

np.save('kminusindexrepo',kminusindex)



# construct mapping between two momenta related by a ky-->-ky rotation
for i in range(npts):
	kx=kxvals[i]
	ky=kyvals[i]
	kxprime=kx
	kyprime=-ky

	if(-kyprime/kxprime>=-1/np.sqrt(3)+10**(-3) and -kyprime/kxprime<1/np.sqrt(3)-10**(-3) and -kxprime>np.sqrt(3)/2*ktheta+10**(-3)):
		kxprime=kxprime+b1x+b2x

		kxprimeprime=kxprime
		kyprimeprime=kyprime
		kxsvalsnew.append(kxprimeprime)
		kysvalsnew.append(kyprimeprime)
	if(-kyprime/kxprime>=-1/np.sqrt(3)+10**(-3) and -kyprime/kxprime<1/np.sqrt(3)+10**(-3) and -kxprime>np.sqrt(3)/2*ktheta+10**(-3)):
		kxprime=kxprime+b1x+b2x

		kxprimeprime=kxprime
		kyprimeprime=kyprime
		kxsvalsnew.append(kxprimeprime)
		kysvalsnew.append(kyprimeprime)


	
	if(-kyprime/kxprime>-1/np.sqrt(3)+10**(-3) and -kyprime/kxprime<1/np.sqrt(3)-10**(-3) and -kxprime<-np.sqrt(3)/2*ktheta+10**(-3) ):
		kxprime=kxprime-b1x-b2x

		kxprimeprime=kxprime
		kyprimeprime=kyprime
		kxsvalsnew.append(kxprimeprime)
		kysvalsnew.append(kyprimeprime)
	if(-kyprime/kxprime>-1/np.sqrt(3)-10**(-3) and -kyprime/kxprime<1/np.sqrt(3)-10**(-3) and -kxprime<-np.sqrt(3)/2*ktheta+10**(-3) ):
		kxprime=kxprime-b1x-b2x

		kxprimeprime=kxprime
		kyprimeprime=kyprime
		kxsvalsnew.append(kxprimeprime)
		kysvalsnew.append(kyprimeprime)




	if(-kxprime/kyprime>=-np.sqrt(3)+10**(-3) and -kxprime<-10**(-3) and kyprime > ktheta-kxprime/np.sqrt(3)+10**(-3)):
		kxprime=kxprime-b1x
		kyprime=kyprime-b1y

		kxprimeprime=kxprime
		kyprimeprime=kyprime

	if(-kxprime/kyprime>=-np.sqrt(3)+10**(-3) and -kxprime<+10**(-3) and kyprime > ktheta-kxprime/np.sqrt(3)+10**(-3)):
		kxprime=kxprime-b1x
		kyprime=kyprime-b1y

		kxprimeprime=kxprime
		kyprimeprime=kyprime



	if(-kxprime/kyprime>=-np.sqrt(3)+10**(-3) and -kxprime>+10**(-3) and kyprime <= -ktheta-kxprime/np.sqrt(3)+10**(-3)):
		kxprime=kxprime+b1x
		kyprime=kyprime+b1y

		kxprimeprime=kxprime
		kyprimeprime=kyprime

	if(-kxprime/kyprime>=-np.sqrt(3)-10**(-3) and -kxprime>+10**(-3) and kyprime <= -ktheta-kxprime/np.sqrt(3)+10**(-3)):
		kxprime=kxprime+b1x
		kyprime=kyprime+b1y

		kxprimeprime=kxprime
		kyprimeprime=kyprime


	

	if(-kxprime/kyprime<np.sqrt(3)-10**(-3) and -kxprime>=10**(-3) and kyprime >= ktheta+kxprime/np.sqrt(3)+10**(-3)):
		kxprime=kxprime+b2x
		kyprime=kyprime+b2y

		kxprimeprime=kxprime
		kyprimeprime=kyprime

	if(-kxprime/kyprime<np.sqrt(3)-10**(-3) and -kxprime>=10**(-3) and kyprime >= ktheta+kxprime/np.sqrt(3)+10**(-3)):
		kxprime=kxprime+b2x
		kyprime=kyprime+b2y

		kxprimeprime=kxprime
		kyprimeprime=kyprime


	if(-kxprime/kyprime<=np.sqrt(3)-10**(-3) and -kxprime<-10**(-3) and kyprime <= -ktheta+kxprime/np.sqrt(3)+10**(-3)):
		kxprime=kxprime-b2x
		kyprime=kyprime-b2y

		kxprimeprime=kxprime
		kyprimeprime=kyprime
	
	if(-kxprime/kyprime<=np.sqrt(3)+10**(-3) and -kxprime<+10**(-3) and kyprime <= -ktheta+kxprime/np.sqrt(3)+10**(-3)):
		kxprime=kxprime-b2x
		kyprime=kyprime-b2y

		kxprimeprime=kxprime
		kyprimeprime=kyprime


	x=np.where(np.round(kxvals,5)==round(kxprime, 5))
	y=np.where(np.round(kyvals,5)==round(kyprime, 5))


	kc2yindex.append(np.intersect1d(x, y)[0])
# REPLACE REFERS TO MY DIRECTORIES, REPLACE WITH OWN PATHS IF USING

np.save('kc2yindexrepo',kc2yindex)

np.save('kxvalsrepo',kxvals)
np.save('kyvalsrepo',kyvals)



# Construct filler matrices for eigenvector matrix
nbands=np.size(qxvalstot)*2
print('NBANDS: ',nbands)
nbandshalf=np.size(qxvalstot)

# This is how many bands we keep, we call it nactive
nbandsminus=nbandshalf-2
nbandsplus=nbandshalf+2
nactive=nbandsplus-nbandsminus
OD=np.zeros([nbands, nactive])
OD2=np.zeros([2*nbands, 2*nactive])


HtotalEigVecs=[]
HtotalEigVals=[]


HtotalEigVecsShell1=[]
HtotalEigVecsShell2=[]
HtotalEigVecsShell3=[]
HtotalEigVecsShell4=[]

print('SHELL 0')

reflection=MZ()

# Here we save wavefunctions for active bands for the first Brillouin zone
for j in range(npts):
	print('pt: ',j)


	Uplus=eighsystem(Hplus(kxvals[j],kyvals[j]))
	Uminus=eighsystem(Hminus(kxvals[j],kyvals[j]))
	ODFull=np.full((nbands,nbands),0)
	Uminustest=np.block([[np.conjugate(eighsystem(Hplus(-kxvals[j],-kyvals[j]))[1]),ODFull],[ODFull,ODFull]])
	Uminusvecs=T0().dot(Uminustest)[nbands:2*nbands,nbandsminus:nbandsplus]

	Htemp=np.block([[Uplus[1][:,nbandsminus:nbandsplus],OD],[OD,Uminusvecs]])
	temp=np.block([[Htemp,OD2],[OD2,Htemp]])

	if j==nptsd-1:
		Uplus=eighsystem(Hplus(kxvals[j],kyvals[j]))
		Uminus=eighsystem(Hminus(kxvals[j],kyvals[j]))
		ODFull=np.full((nbands,nbands),0)

		Uminustest=np.block([[np.conjugate(eighsystem(Hplus(-kxvals[j],-kyvals[j]))[1]),ODFull],[ODFull,ODFull]])
		Uminusvecs=T0().dot(Uminustest)[nbands:2*nbands,nbandsminus:nbandsplus]




		tempplus=Uplus[1][:,nbandsminus:nbandsplus]

		if np.abs(D)<10**(-3)*wab12:
			# We do a qr decomposition of the default eigenvectors after action of C_{2zT} to get an orthornomal set that respects C2zT

			tempplus[:,[0,1]] = tempplus[:,[1,0]]
			tempplus[:,[2,3]] = tempplus[:,[3,2]]
			tempplus[:,0]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempplus[:,0]))+tempplus[:,0])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempplus[:,0]))+tempplus[:,0])
			tempplus[:,1]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempplus[:,1]))+tempplus[:,1])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempplus[:,1]))+tempplus[:,1])
			tempplus[:,2]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempplus[:,2]))+tempplus[:,2])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempplus[:,2]))+tempplus[:,2])
			tempplus[:,3]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempplus[:,3]))+tempplus[:,3])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempplus[:,3]))+tempplus[:,3])
		
		q,r=np.linalg.qr(tempplus)
		tempplus=q


		Htemp=np.block([[tempplus,OD],[OD,Uminusvecs]])


		


	if j==0:
		Uplus=eighsystem(Hplus(kxvals[j],kyvals[j]))
		Uminus=eighsystem(Hminus(kxvals[j],kyvals[j]))
		ODFull=np.full((nbands,nbands),0)

		Uminustest=np.block([[np.conjugate(eighsystem(Hplus(-kxvals[j],-kyvals[j]))[1]),ODFull],[ODFull,ODFull]])
		Uminusvecs=T0().dot(Uminustest)[nbands:2*nbands,nbandsminus:nbandsplus]




		tempminus=Uminusvecs

		if np.abs(D)<10**(-3)*wab12:
			tempminus[:,[0,1]] = tempminus[:,[1,0]]
			tempminus[:,[3,2]] = tempminus[:,[2,3]]
			tempminus[:,0]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempminus[:,0]))+tempminus[:,0])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempminus[:,0]))+tempminus[:,0])
			tempminus[:,1]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempminus[:,1]))+tempminus[:,1])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempminus[:,1]))+tempminus[:,1])
			tempminus[:,2]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempminus[:,2]))+tempminus[:,2])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempminus[:,2]))+tempminus[:,2])
			tempminus[:,3]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempminus[:,3]))+tempminus[:,3])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempminus[:,3]))+tempminus[:,3])
		q,r=np.linalg.qr(tempminus)
		tempminus=q
		

		Htemp=np.block([[Uplus[1][:,nbandsminus:nbandsplus],OD],[OD,tempminus]])



	kin=[]

	for i in range(nactive):

		kin.append(Uplus[0][nbandsminus+i])

	for i in range(nactive):

		kin.append(Uminus[0][nbandsminus+i])

	for i in range(nactive):

		kin.append(Uplus[0][nbandsminus+i])

	for i in range(nactive):

		kin.append(Uminus[0][nbandsminus+i])


	if np.abs(kxvals[j])<10**(-4)*ktheta and np.abs(kyvals[j])<10**(-4)*ktheta:
		print('correcting')

		vecsplusgood=Uplus[1][:,nbandsminus+1:nbandsplus-1]
		vecsplushigh=Uplus[1][:,nbandsplus-1:nbandsplus+1]
		vecspluslow=Uplus[1][:,nbandsminus-1:nbandsminus+1]

		hcorrectedplus=np.hstack((vecspluslow,vecsplusgood,vecsplushigh))[:,1:nactive+1]

		hcorrectedplus[:,0]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(hcorrectedplus[:,0]))+hcorrectedplus[:,0])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(hcorrectedplus[:,0]))+hcorrectedplus[:,0])
		hcorrectedplus[:,1]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(hcorrectedplus[:,1]))+hcorrectedplus[:,1])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(hcorrectedplus[:,1]))+hcorrectedplus[:,1])
		hcorrectedplus[:,2]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(hcorrectedplus[:,2]))+hcorrectedplus[:,2])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(hcorrectedplus[:,2]))+hcorrectedplus[:,2])
		hcorrectedplus[:,3]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(hcorrectedplus[:,3]))+hcorrectedplus[:,3])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(hcorrectedplus[:,3]))+hcorrectedplus[:,3])

		q,r=np.linalg.qr(hcorrectedplus)
		hcorrectedplus=q


		hcorrectedminus=(T0().dot(np.block([[np.conjugate(hcorrectedplus),OD],[OD,OD]])))[nbands:2*nbands,0:nactive]




		Htemp=np.block([[hcorrectedplus,OD],[OD,hcorrectedminus]])



	Htemp=Htemp.dot(np.sqrt(np.diag(np.diagonal(np.conjugate(np.transpose(Htemp)).dot(C2zT()).dot(np.conjugate(Htemp))))))
	



	Htemp[:,2] *= np.sign(np.imag(np.conjugate(np.transpose(Htemp)).dot(Chiral()).dot(Htemp)[2][1]))
	Htemp[:,6] *= -np.sign(np.imag(np.conjugate(np.transpose(Htemp)).dot(Chiral()).dot(Htemp)[6][5]))



	Htemp[:,5] *= -np.sign(np.real(np.conjugate(np.transpose(Htemp)).dot(PHC2z()).dot(Htemp))[2][5])
	Htemp[:,6] *= np.sign(np.real(np.conjugate(np.transpose(Htemp)).dot(PHC2z()).dot(Htemp))[1][6])

	temp=np.block([[Htemp,OD2],[OD2,Htemp]])



	HtotalEigVecs.append(np.block([[Htemp,OD2],[OD2,Htemp]]))
	HtotalEigVals.append(kin)
# REPLACE REFERS TO MY DIRECTORIES, REPLACE WITH OWN PATHS IF USING

np.savez_compressed('/n/holylfs/LABS/sachdev_lab/vectors2/HtotalEigVecsrepo_%d'%(int(sys.argv[1]),),HtotalEigVecs)
np.savez_compressed('/n/holylfs/LABS/sachdev_lab/vectors2/HtotalEigValsrepo_%d'%(int(sys.argv[1]),),HtotalEigVals)
print('SHELL 1')

# Here we save wavefunctions for active bands in the first shell of Moire unit cells
for i in range(np.size(shell1momentax)):
	print('shell: ',i)

	HtotalEigVecsShell1x=[]
	for j in range(npts):
		print('pt: ',j)


	
		ODFull=np.full((nbands,nbands),0)
	

		Uplus=eighsystem(Hplus(kxvals[j]+shell1momentax[i],kyvals[j]+shell1momentay[i]))
		Uminus=eighsystem(Hminus(kxvals[j]+shell1momentax[i],kyvals[j]+shell1momentay[i]))
		Uminustest=np.block([[np.conjugate(eighsystem(Hplus(-kxvals[j]-shell1momentax[i],-kyvals[j]-shell1momentay[i]))[1]),ODFull],[ODFull,ODFull]])
		Uminusvecs=T0().dot(Uminustest)[nbands:2*nbands,nbandsminus:nbandsplus]


		unew=np.full(nbands,0.+0.j)
		hplusnew=np.full((nbands,nactive),0.+0.j)
		hminusnew=np.full((nbands,nactive),0.+0.j)




		Htemp=np.block([[Uplus[1][:,nbandsminus:nbandsplus],OD],[OD,Uminusvecs]])

		if j==nptsd-1:
			Uplus=eighsystem(Hplus(kxvals[j]+shell1momentax[i],kyvals[j]+shell1momentay[i]))
			Uminus=eighsystem(Hminus(kxvals[j]+shell1momentax[i],kyvals[j]+shell1momentay[i]))
			ODFull=np.full((nbands,nbands),0)

			Uminustest=np.block([[np.conjugate(eighsystem(Hplus(-kxvals[j]-shell1momentax[i],-kyvals[j]-shell1momentay[i]))[1]),ODFull],[ODFull,ODFull]])
			Uminusvecs=T0().dot(Uminustest)[nbands:2*nbands,nbandsminus:nbandsplus]




			tempplus=eighsystem(Hplus(kxvals[j]+shell1momentax[i],kyvals[j]+shell1momentay[i]))[1][:,nbandsminus:nbandsplus]

			# print(tempminus)
			if np.abs(D)<10**(-3)*wab12:
				tempplus[:,[0,1]] = tempplus[:,[1,0]]
				tempplus[:,[2,3]] = tempplus[:,[3,2]]
				tempplus[:,0]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempplus[:,0]))+tempplus[:,0])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempplus[:,0]))+tempplus[:,0])
				tempplus[:,1]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempplus[:,1]))+tempplus[:,1])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempplus[:,1]))+tempplus[:,1])
				tempplus[:,2]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempplus[:,2]))+tempplus[:,2])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempplus[:,2]))+tempplus[:,2])
				tempplus[:,3]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempplus[:,3]))+tempplus[:,3])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempplus[:,3]))+tempplus[:,3])
				q,r=np.linalg.qr(tempplus)
				tempplus=q


			Htemp=np.block([[tempplus,OD],[OD,Uminusvecs]])

			


		if j==0:
			Uplus=eighsystem(Hplus(kxvals[j]+shell1momentax[i],kyvals[j]+shell1momentay[i]))
			Uminus=eighsystem(Hminus(kxvals[j]+shell1momentax[i],kyvals[j]+shell1momentay[i]))
			ODFull=np.full((nbands,nbands),0)

			Uminustest=np.block([[np.conjugate(eighsystem(Hplus(-kxvals[j]-shell1momentax[i],-kyvals[j]-shell1momentay[i]))[1]),ODFull],[ODFull,ODFull]])
			Uminusvecs=T0().dot(Uminustest)[nbands:2*nbands,nbandsminus:nbandsplus]




			tempminus=Uminusvecs

			# print(tempminus)
			if np.abs(D)<10**(-3)*wab12:
				tempminus[:,[0,1]] = tempminus[:,[1,0]]
				tempminus[:,[2,3]] = tempminus[:,[3,2]]
				tempminus[:,0]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempminus[:,0]))+tempminus[:,0])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempminus[:,0]))+tempminus[:,0])
				tempminus[:,1]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempminus[:,1]))+tempminus[:,1])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempminus[:,1]))+tempminus[:,1])
				tempminus[:,2]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempminus[:,2]))+tempminus[:,2])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempminus[:,2]))+tempminus[:,2])
				tempminus[:,3]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempminus[:,3]))+tempminus[:,3])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempminus[:,3]))+tempminus[:,3])
				q,r=np.linalg.qr(tempminus)
				tempminus=q




			Htemp=np.block([[Uplus[1][:,nbandsminus:nbandsplus],OD],[OD,tempminus]])




		if np.abs(kxvals[j])<10**(-4)*ktheta and np.abs(kyvals[j])<10**(-4)*ktheta:
			print('correcting')

			vecsplusgood=Uplus[1][:,nbandsminus+1:nbandsplus-1]
			vecsplushigh=Uplus[1][:,nbandsplus-1:nbandsplus+1]
			vecspluslow=Uplus[1][:,nbandsminus-1:nbandsminus+1]

			hcorrectedplus=np.hstack((vecspluslow,vecsplusgood,vecsplushigh))[:,1:nactive+1]

			hcorrectedplus[:,0]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(hcorrectedplus[:,0]))+hcorrectedplus[:,0])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(hcorrectedplus[:,0]))+hcorrectedplus[:,0])
			hcorrectedplus[:,1]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(hcorrectedplus[:,1]))+hcorrectedplus[:,1])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(hcorrectedplus[:,1]))+hcorrectedplus[:,1])
			hcorrectedplus[:,2]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(hcorrectedplus[:,2]))+hcorrectedplus[:,2])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(hcorrectedplus[:,2]))+hcorrectedplus[:,2])
			hcorrectedplus[:,3]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(hcorrectedplus[:,3]))+hcorrectedplus[:,3])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(hcorrectedplus[:,3]))+hcorrectedplus[:,3])

			q,r=np.linalg.qr(hcorrectedplus)
			hcorrectedplus=q


			hcorrectedminus=(T0().dot(np.block([[np.conjugate(hcorrectedplus),OD],[OD,OD]])))[nbands:2*nbands,0:nactive]




			Htemp=np.block([[hcorrectedplus,OD],[OD,hcorrectedminus]])





		Htemp=Htemp.dot(np.sqrt(np.diag(np.diagonal(np.conjugate(np.transpose(Htemp)).dot(C2zT()).dot(np.conjugate(Htemp))))))

		Htemp[:,2] *= np.sign(np.imag(np.conjugate(np.transpose(Htemp)).dot(Chiral()).dot(Htemp)[2][1]))
		Htemp[:,6] *= -np.sign(np.imag(np.conjugate(np.transpose(Htemp)).dot(Chiral()).dot(Htemp)[6][5]))



		Htemp[:,5] *= -np.sign(np.real(np.conjugate(np.transpose(Htemp)).dot(PHC2z()).dot(Htemp)))[2][5]
		Htemp[:,6] *= np.sign(np.real(np.conjugate(np.transpose(Htemp)).dot(PHC2z()).dot(Htemp)))[1][6]


		HtotalEigVecsShell1x.append(np.block([[Htemp,OD2],[OD2,Htemp]]))
	HtotalEigVecsShell1.append(HtotalEigVecsShell1x)


HtotalEigVecsShell1newnew=np.full((np.size(shell1momentax),npts,4*nbands,4*nactive),0.0j)
HtotalEigVecscopy=HtotalEigVecs
HtotalEigVecsShell1copy=HtotalEigVecsShell1

print('SHELL 1')
# Fix sign of wavefunction in different moire unit cells to enforce periodic gauge. For a model constructed with 3 shells,we find the waveufnctions are only close enough to periodic to fix the sign in the first shell of moire unit cells. However the gauge fixing condition will effect the energies only at the level of the Hartree term so we will only use the first shell where the gauge is fixed for the Hartree contribution (else the hermiticity constraints will begin to break down and the sum over form factor components which should be odd will not cancel). This sectin is just a check as printed in the next loop. We will fix the gauge for real after we rotate the wavefunctions later.
for i in range(np.size(shell1momentax)):
	print('shell: ',i)
	HtotalEigVecsShell1x=[]
	for j in range(npts):

		countminus=0
		countplus=0
		countcheck=0

		for k in range(np.size(qvalslayer)):

			checkplusx=np.where(np.round(qxvalstot/ktheta,3)==np.round((qxvalstot[k]-shell1momentax[i])/ktheta,3))
			checkplusy=np.where(np.round(qyvalstot/ktheta,3)==np.round((qyvalstot[k]-shell1momentay[i])/ktheta,3))

			checkpluslayer=np.where(np.round(qvalslayer,3)==np.round(qvalslayer[k],3))
			intersectplus=np.intersect1d(checkpluslayer,np.intersect1d(checkplusx, checkplusy))


			checkminusx=np.where(np.round(qxvalstotminus/ktheta,3)==np.round((qxvalstotminus[k]-shell1momentax[i])/ktheta,3))
			checkminusy=np.where(np.round(qyvalstotminus/ktheta,3)==np.round((qyvalstotminus[k]-shell1momentay[i])/ktheta,3))

			checkminuslayer=np.where(np.round(qvalslayer,3)==np.round(qvalslayer[k],3))
			intersectminus=np.intersect1d(checkminuslayer,np.intersect1d(checkminusx, checkminusy))
			if np.size(intersectminus)>0  and countminus==0:
				if np.abs(HtotalEigVecscopy[j][nbands+2*intersectminus[0],4])<10**(-10):
					continue
				countminus=countminus+1



				phaserowv2up=HtotalEigVecscopy[j][nbands+2*intersectminus[0],nactive:2*nactive]/HtotalEigVecsShell1copy[i][j][nbands+2*k,nactive:2*nactive]*np.abs(HtotalEigVecsShell1copy[i][j][nbands+2*k,nactive:2*nactive])/np.abs(HtotalEigVecscopy[j][nbands+2*intersectminus[0],nactive:2*nactive])

				phaserowv2down=HtotalEigVecscopy[j][3*nbands+2*intersectminus[0],3*nactive:4*nactive]/HtotalEigVecsShell1copy[i][j][3*nbands+2*k,3*nactive:4*nactive]*np.abs(HtotalEigVecsShell1copy[i][j][3*nbands+2*k,3*nactive:4*nactive])/np.abs(HtotalEigVecscopy[j][3*nbands+2*intersectminus[0],3*nactive:4*nactive])



			if np.size(intersectplus)>0  and countplus==0:
				if np.abs(HtotalEigVecscopy[j][2*intersectplus[0],0])<10**(-10) :
					continue
				countplus=countplus+1

				phaserowv1up=HtotalEigVecscopy[j][2*intersectplus[0],0:nactive]/HtotalEigVecsShell1copy[i][j][2*k,0:nactive]*np.abs(HtotalEigVecsShell1copy[i][j][2*k,0:nactive])/np.abs(HtotalEigVecscopy[j][2*intersectplus[0],0:nactive])
				
				phaserowv1down=HtotalEigVecscopy[j][2*nbands+2*intersectplus[0],2*nactive:3*nactive]/HtotalEigVecsShell1copy[i][j][2*nbands+2*k,2*nactive:3*nactive]*np.abs(HtotalEigVecsShell1copy[i][j][2*nbands+2*k,2*nactive:3*nactive])/np.abs(HtotalEigVecscopy[j][2*nbands+2*intersectplus[0],2*nactive:3*nactive])

				
			if countplus>0 and countminus>0:
				

				phaserow=np.sign(np.real((np.hstack((phaserowv1up,phaserowv2up,phaserowv1down,phaserowv2down)))))

				countcheck=countcheck+1




				phase=np.diag(phaserow)
				


				HtotalEigVecsShell1newnew[i][j]=HtotalEigVecsShell1copy[i][j].dot(phase)



# REPLACE REFERS TO MY DIRECTORIES, REPLACE WITH OWN PATHS IF USING

np.savez_compressed('/n/holylfs/LABS/sachdev_lab/vectors2/HtotalEigVecsShell1repo_%d'%(int(sys.argv[1]),),HtotalEigVecsShell1)
print('SHELL 2')
# Here we save wavefunctions for active bands in the first shell of Moire unit cells
for i in range(np.size(shell2momentax)):
	print('shell: ',i)

	HtotalEigVecsShell2x=[]
	for j in range(npts):
		print('pt: ',j)


	
		ODFull=np.full((nbands,nbands),0)
	

		Uplus=eighsystem(Hplus(kxvals[j]+shell2momentax[i],kyvals[j]+shell2momentay[i]))
		Uminus=eighsystem(Hminus(kxvals[j]+shell2momentax[i],kyvals[j]+shell2momentay[i]))
		Uminustest=np.block([[np.conjugate(eighsystem(Hplus(-kxvals[j]-shell2momentax[i],-kyvals[j]-shell2momentay[i]))[1]),ODFull],[ODFull,ODFull]])
		Uminusvecs=T0().dot(Uminustest)[nbands:2*nbands,nbandsminus:nbandsplus]


		unew=np.full(nbands,0.+0.j)
		hplusnew=np.full((nbands,nactive),0.+0.j)
		hminusnew=np.full((nbands,nactive),0.+0.j)


		Uplusvecs=HtotalEigVecs[j][0:nbands,0:nactive]
		Uminusvecs=HtotalEigVecs[j][nbands:2*nbands,nactive:2*nactive]



		Htemp=np.block([[Uplus[1][:,nbandsminus:nbandsplus],OD],[OD,Uminusvecs]])

		if j==nptsd-1:
			Uplus=eighsystem(Hplus(kxvals[j]+shell2momentax[i],kyvals[j]+shell2momentay[i]))
			Uminus=eighsystem(Hminus(kxvals[j]+shell2momentax[i],kyvals[j]+shell2momentay[i]))
			ODFull=np.full((nbands,nbands),0)

			Uminustest=np.block([[np.conjugate(eighsystem(Hplus(-kxvals[j]-shell2momentax[i],-kyvals[j]-shell2momentay[i]))[1]),ODFull],[ODFull,ODFull]])
			Uminusvecs=T0().dot(Uminustest)[nbands:2*nbands,nbandsminus:nbandsplus]




			tempplus=eighsystem(Hplus(kxvals[j]+shell2momentax[i],kyvals[j]+shell2momentay[i]))[1][:,nbandsminus:nbandsplus]

			if np.abs(D)<10**(-3)*wab12:
				tempplus[:,[0,1]] = tempplus[:,[1,0]]
				tempplus[:,[2,3]] = tempplus[:,[3,2]]
				tempplus[:,0]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempplus[:,0]))+tempplus[:,0])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempplus[:,0]))+tempplus[:,0])
				tempplus[:,1]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempplus[:,1]))+tempplus[:,1])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempplus[:,1]))+tempplus[:,1])
				tempplus[:,2]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempplus[:,2]))+tempplus[:,2])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempplus[:,2]))+tempplus[:,2])
				tempplus[:,3]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempplus[:,3]))+tempplus[:,3])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempplus[:,3]))+tempplus[:,3])

			q,r=np.linalg.qr(tempplus)
			tempplus=q

			Htemp=np.block([[tempplus,OD],[OD,Uminusvecs]])

			


		if j==0:
			Uplus=eighsystem(Hplus(kxvals[j]+shell2momentax[i],kyvals[j]+shell2momentay[i]))
			Uminus=eighsystem(Hminus(kxvals[j]+shell2momentax[i],kyvals[j]+shell2momentay[i]))
			ODFull=np.full((nbands,nbands),0)

			Uminustest=np.block([[np.conjugate(eighsystem(Hplus(-kxvals[j]-shell2momentax[i],-kyvals[j]-shell2momentay[i]))[1]),ODFull],[ODFull,ODFull]])
			Uminusvecs=T0().dot(Uminustest)[nbands:2*nbands,nbandsminus:nbandsplus]




			tempminus=Uminusvecs

			if np.abs(D)<10**(-3)*wab12:
				tempminus[:,[0,1]] = tempminus[:,[1,0]]
				tempminus[:,[2,3]] = tempminus[:,[3,2]]
				tempminus[:,0]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempminus[:,0]))+tempminus[:,0])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempminus[:,0]))+tempminus[:,0])
				tempminus[:,1]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempminus[:,1]))+tempminus[:,1])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempminus[:,1]))+tempminus[:,1])
				tempminus[:,2]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempminus[:,2]))+tempminus[:,2])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempminus[:,2]))+tempminus[:,2])
				tempminus[:,3]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempminus[:,3]))+tempminus[:,3])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(tempminus[:,3]))+tempminus[:,3])
			q,r=np.linalg.qr(tempminus)
			tempminus=q



			Htemp=np.block([[Uplus[1][:,nbandsminus:nbandsplus],OD],[OD,tempminus]])


		if np.abs(kxvals[j])<10**(-4)*ktheta and np.abs(kyvals[j])<10**(-4)*ktheta:
			print('correcting')

			vecsplusgood=Uplus[1][:,nbandsminus+1:nbandsplus-1]
			vecsplushigh=Uplus[1][:,nbandsplus-1:nbandsplus+1]
			vecspluslow=Uplus[1][:,nbandsminus-1:nbandsminus+1]

			hcorrectedplus=np.hstack((vecspluslow,vecsplusgood,vecsplushigh))[:,1:nactive+1]

			hcorrectedplus[:,0]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(hcorrectedplus[:,0]))+hcorrectedplus[:,0])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(hcorrectedplus[:,0]))+hcorrectedplus[:,0])
			hcorrectedplus[:,1]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(hcorrectedplus[:,1]))+hcorrectedplus[:,1])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(hcorrectedplus[:,1]))+hcorrectedplus[:,1])
			hcorrectedplus[:,2]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(hcorrectedplus[:,2]))+hcorrectedplus[:,2])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(hcorrectedplus[:,2]))+hcorrectedplus[:,2])
			hcorrectedplus[:,3]=(C2zT()[0:nbands,0:nbands].dot(np.conjugate(hcorrectedplus[:,3]))+hcorrectedplus[:,3])/np.linalg.norm(C2zT()[0:nbands,0:nbands].dot(np.conjugate(hcorrectedplus[:,3]))+hcorrectedplus[:,3])

			q,r=np.linalg.qr(hcorrectedplus)
			hcorrectedplus=q


			hcorrectedminus=(T0().dot(np.block([[np.conjugate(hcorrectedplus),OD],[OD,OD]])))[nbands:2*nbands,0:nactive]





			Htemp=np.block([[hcorrectedplus,OD],[OD,hcorrectedminus]])




		Htemp=Htemp.dot(np.sqrt(np.diag(np.diagonal(np.conjugate(np.transpose(Htemp)).dot(C2zT()).dot(np.conjugate(Htemp))))))

		Htemp[:,2] *= np.sign(np.imag(np.conjugate(np.transpose(Htemp)).dot(Chiral()).dot(Htemp)[2][1]))
		Htemp[:,6] *= -np.sign(np.imag(np.conjugate(np.transpose(Htemp)).dot(Chiral()).dot(Htemp)[6][5]))



		Htemp[:,5] *= -np.sign(np.real(np.conjugate(np.transpose(Htemp)).dot(PHC2z()).dot(Htemp)))[2][5]
		Htemp[:,6] *= np.sign(np.real(np.conjugate(np.transpose(Htemp)).dot(PHC2z()).dot(Htemp)))[1][6]





		HtotalEigVecsShell2x.append(np.block([[Htemp,OD2],[OD2,Htemp]]))
	HtotalEigVecsShell2.append(HtotalEigVecsShell2x)
# REPLACE REFERS TO MY DIRECTORIES, REPLACE WITH OWN PATHS IF USING

np.savez_compressed('/n/holylfs/LABS/sachdev_lab/vectors2/HtotalEigVecsShell2repo_%d'%(int(sys.argv[1]),),HtotalEigVecsShell2)

HtotalEigVecscopy=HtotalEigVecs
HtotalEigVecsShell1copy=HtotalEigVecsShell1
HtotalEigVecsShell2copy=HtotalEigVecsShell2
HtotalEigVecsShell3copy=HtotalEigVecsShell3
HtotalEigVecsShell4copy=HtotalEigVecsShell4

HtotalEigVecs=[]
HtotalEigVecsShell1=[]
HtotalEigVecsShell2=[]
HtotalEigVecsShell3=[]
HtotalEigVecsShell4=[]


print('SHELL 0')

#Here we use C3 and C2z to make sure these symmetries of the eigenevectors are exact. While they are to good approximation without rotating the eigenvectors, we need these symmetries respected to errors of 10^(-10) in order to remain symmetries under HF iterationns, reqiring rotations
for j in range(npts):
	print('pt: ',j)

	blockzer=np.full((nbands,nactive),0.0+0.0j)
	blockzerbig=np.full((nbands,nactive),0.0+0.0j)

	if np.abs(kxvals[j]/ktheta)<10**(-3) and np.abs(kyvals[j]/ktheta)<10**(-3) :

		hplus=HtotalEigVecscopy[j][0:nbands,0:nactive]
		kxprime,kyprime=C2zrotate(kxvals[j],kyvals[j])
		kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
		kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
		indexfind=np.intersect1d(kxfind,kyfind)

		if np.size(indexfind)==0:
			for k in range(np.size(shell1momentax)):
				kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
				kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
				indexfind=np.intersect1d(kxfind,kyfind)
				if np.size(indexfind)>0:
					indexnew=indexfind[0]
					kxprime=kxvals[indexnew]+shell1momentax[k]
					kyprime=kyvals[indexnew]+shell1momentay[k]
					U=HtotalEigVecsShell1copy[k][indexnew][0:nbands,0:nactive]
		else:
			indexnew=indexfind[0]
			kxprime=kxvals[indexnew]
			kyprime=kyvals[indexnew]
			U=HtotalEigVecscopy[indexnew][0:nbands,0:nactive]
		
		hminus=C2zfull().dot(np.block([[U,blockzer],[blockzer,blockzer]]))[nbands:2*nbands,0:nactive]




		Htemp=np.block([[hplus,blockzerbig],[blockzerbig,hminus]])


	


	elif kxvals[j]/ktheta<0-10**(-3) and kxvals[j]/ktheta*q3y/q3x-kyvals[j]/ktheta<10**(-3) :

		hplus=HtotalEigVecscopy[j][0:nbands,0:nactive]
		kxprime,kyprime=C2zrotate(kxvals[j],kyvals[j])
		kxprime,kyprime=C3rotate(kxprime,kyprime)
		kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
		kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
		indexfind=np.intersect1d(kxfind,kyfind)

		if np.size(indexfind)==0:
			for k in range(np.size(shell1momentax)):
				kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
				kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
				indexfind=np.intersect1d(kxfind,kyfind)
				if np.size(indexfind)>0:
					indexnew=indexfind[0]
					kxprime=kxvals[indexnew]+shell1momentax[k]
					kyprime=kyvals[indexnew]+shell1momentay[k]
					U=HtotalEigVecsShell1copy[k][indexnew][0:nbands,0:nactive]
		else:
			indexnew=indexfind[0]
			kxprime=kxvals[indexnew]
			kyprime=kyvals[indexnew]
			U=HtotalEigVecscopy[indexnew][0:nbands,0:nactive]
		
		hminus=(C3full()).dot(C3full()).dot(C2zfull()).dot(np.block([[U,blockzer],[blockzer,blockzer]]))[nbands:2*nbands,0:nactive]


		Htemp=np.block([[hplus,blockzerbig],[blockzerbig,hminus]])


	elif kxvals[j]/ktheta>0-10**(-3) and -kxvals[j]/ktheta*q3y/q3x-kyvals[j]/ktheta<-10**(-3):

		hplus=HtotalEigVecscopy[j][0:nbands,0:nactive]
		kxprime,kyprime=C2zrotate(kxvals[j],kyvals[j])
		kxprime,kyprime=C3rotate(kxprime,kyprime)
		kxprime,kyprime=C3rotate(kxprime,kyprime)
		kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
		kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
		indexfind=np.intersect1d(kxfind,kyfind)

		if np.size(indexfind)==0:
			for k in range(np.size(shell1momentax)):
				kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
				kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
				indexfind=np.intersect1d(kxfind,kyfind)
				if np.size(indexfind)>0:
					indexnew=indexfind[0]
					kxprime=kxvals[indexnew]+shell1momentax[k]
					kyprime=kyvals[indexnew]+shell1momentay[k]
					U=HtotalEigVecsShell1copy[k][indexnew][0:nbands,0:nactive]
		else:
			indexnew=indexfind[0]
			kxprime=kxvals[indexnew]
			kyprime=kyvals[indexnew]
			U=HtotalEigVecscopy[indexnew][0:nbands,0:nactive]
		
		hminus=(C3full()).dot(C2zfull()).dot(np.block([[U,blockzer],[blockzer,blockzer]]))[nbands:2*nbands,0:nactive]


		Htemp=np.block([[hplus,blockzerbig],[blockzerbig,hminus]])




	elif kxvals[j]/ktheta>0+10**(-3) and -kxvals[j]/ktheta*q3y/q3x-kyvals[j]/ktheta>-10**(-3) and kxvals[j]/ktheta*q3y/q3x-kyvals[j]/ktheta<-10**(-3):
		kxprime,kyprime=C3rotate(kxvals[j],kyvals[j])

		kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
		kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
		indexfind=np.intersect1d(kxfind,kyfind)[0]

		# print(kxprime/ktheta,kyprime/ktheta)
		kxprime=kxvals[indexfind]
		kyprime=kyvals[indexfind]

		U=HtotalEigVecscopy[indexfind][0:2*nbands,0:2*nactive]
		hplus=C3full().dot(C3full()).dot(U)[0:nbands,0:nactive]



		kxprime,kyprime=C2zrotate(kxvals[j],kyvals[j])
		kxprime,kyprime=C3rotate(kxprime,kyprime)
		kxprime,kyprime=C3rotate(kxprime,kyprime)
		kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
		kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
		indexfind=np.intersect1d(kxfind,kyfind)

		if np.size(indexfind)==0:
			for k in range(np.size(shell1momentax)):
				kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
				kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
				indexfind=np.intersect1d(kxfind,kyfind)
				if np.size(indexfind)>0:
					indexnew=indexfind[0]
					kxprime=kxvals[indexnew]+shell1momentax[k]
					kyprime=kyvals[indexnew]+shell1momentay[k]
					U=HtotalEigVecsShell1copy[k][indexnew][0:nbands,0:nactive]
		else:
			indexnew=indexfind[0]
			kxprime=kxvals[indexnew]
			kyprime=kyvals[indexnew]
			U=HtotalEigVecscopy[indexnew][0:nbands,0:nactive]
		
		hminus=(C3full()).dot(C2zfull()).dot(np.block([[U,blockzer],[blockzer,blockzer]]))[nbands:2*nbands,0:nactive]


		Htemp=np.block([[hplus,blockzerbig],[blockzerbig,hminus]])


	elif kxvals[j]/ktheta>0+10**(-3) and kxvals[j]/ktheta*q3y/q3x-kyvals[j]/ktheta>-10**(-3):
		kxprime,kyprime=C3rotate(kxvals[j],kyvals[j])

		kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
		kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
		indexfind=np.intersect1d(kxfind,kyfind)[0]

		kxprime=kxvals[indexfind]
		kyprime=kyvals[indexfind]

		U=HtotalEigVecscopy[indexfind][0:2*nbands,0:2*nactive]
		hplus=C3full().dot(C3full()).dot(U)[0:nbands,0:nactive]

		kxprime,kyprime=C2zrotate(kxvals[j],kyvals[j])
		kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
		kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
		indexfind=np.intersect1d(kxfind,kyfind)

		if np.size(indexfind)==0:
			for k in range(np.size(shell1momentax)):
				kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
				kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
				indexfind=np.intersect1d(kxfind,kyfind)
				if np.size(indexfind)>0:
					indexnew=indexfind[0]
					kxprime=kxvals[indexnew]+shell1momentax[k]
					kyprime=kyvals[indexnew]+shell1momentay[k]
					U=HtotalEigVecsShell1copy[k][indexnew][0:nbands,0:nactive]
		else:
			indexnew=indexfind[0]
			kxprime=kxvals[indexnew]
			kyprime=kyvals[indexnew]
			U=HtotalEigVecscopy[indexnew][0:nbands,0:nactive]
		
		hminus=C2zfull().dot(np.block([[U,blockzer],[blockzer,blockzer]]))[nbands:2*nbands,0:nactive]


		Htemp=np.block([[hplus,blockzerbig],[blockzerbig,hminus]])


	elif kxvals[j]/ktheta<0+10**(-3) and -kxvals[j]/ktheta*q3y/q3x-kyvals[j]/ktheta>10**(-3):
		kxprime,kyprime=C3rotate(kxvals[j],kyvals[j])
		kxprime,kyprime=C3rotate(kxprime,kyprime)

		kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
		kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
		indexfind=np.intersect1d(kxfind,kyfind)[0]

		kxprime=kxvals[indexfind]
		kyprime=kyvals[indexfind]

		U=HtotalEigVecscopy[indexfind][0:2*nbands,0:2*nactive]
		hplus=C3full().dot(U)[0:nbands,0:nactive]
		kxprime,kyprime=C2zrotate(kxvals[j],kyvals[j])
		kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
		kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
		indexfind=np.intersect1d(kxfind,kyfind)

		if np.size(indexfind)==0:
			for k in range(np.size(shell1momentax)):
				kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
				kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
				indexfind=np.intersect1d(kxfind,kyfind)
				if np.size(indexfind)>0:
					indexnew=indexfind[0]
					kxprime=kxvals[indexnew]+shell1momentax[k]
					kyprime=kyvals[indexnew]+shell1momentay[k]
					U=HtotalEigVecsShell1copy[k][indexnew][0:nbands,0:nactive]
		else:
			indexnew=indexfind[0]
			kxprime=kxvals[indexnew]
			kyprime=kyvals[indexnew]
			U=HtotalEigVecscopy[indexnew][0:nbands,0:nactive]
		
		hminus=C2zfull().dot(np.block([[U,blockzer],[blockzer,blockzer]]))[nbands:2*nbands,0:nactive]


		Htemp=np.block([[hplus,blockzerbig],[blockzerbig,hminus]])


	elif kxvals[j]/ktheta<0-10**(-3) and -kxvals[j]/ktheta*q3y/q3x-kyvals[j]/ktheta<10**(-3) and kxvals[j]/ktheta*q3y/q3x-kyvals[j]/ktheta>10**(-3):
		kxprime,kyprime=C3rotate(kxvals[j],kyvals[j])
		kxprime,kyprime=C3rotate(kxprime,kyprime)

		kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
		kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
		indexfind=np.intersect1d(kxfind,kyfind)[0]

		kxprime=kxvals[indexfind]
		kyprime=kyvals[indexfind]


		U=HtotalEigVecscopy[indexfind][0:2*nbands,0:2*nactive]
		hplus=C3full().dot(U)[0:nbands,0:nactive]
		kxprime,kyprime=C2zrotate(kxvals[j],kyvals[j])
		kxprime,kyprime=C3rotate(kxprime,kyprime)
		kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
		kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
		indexfind=np.intersect1d(kxfind,kyfind)

		if np.size(indexfind)==0:
			for k in range(np.size(shell1momentax)):
				kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
				kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
				indexfind=np.intersect1d(kxfind,kyfind)
				if np.size(indexfind)>0:
					indexnew=indexfind[0]
					kxprime=kxvals[indexnew]+shell1momentax[k]
					kyprime=kyvals[indexnew]+shell1momentay[k]
					U=HtotalEigVecsShell1copy[k][indexnew][0:nbands,0:nactive]
		else:
			indexnew=indexfind[0]
			kxprime=kxvals[indexnew]
			kyprime=kyvals[indexnew]
			U=HtotalEigVecscopy[indexnew][0:nbands,0:nactive]
		
		hminus=(C3full()).dot(C3full()).dot(C2zfull()).dot(np.block([[U,blockzer],[blockzer,blockzer]]))[nbands:2*nbands,0:nactive]


		Htemp=np.block([[hplus,blockzerbig],[blockzerbig,hminus]])



	Htemp[:,4] *= -np.sign(np.real(np.conjugate(np.transpose(Htemp)).dot(PHC2z()).dot(Htemp)[2][5]))
	Htemp[:,7] *= np.sign(np.real(np.conjugate(np.transpose(Htemp)).dot(PHC2z()).dot(Htemp)[1][6]))

	Htemp[:,5] *= -np.sign(np.real(np.conjugate(np.transpose(Htemp)).dot(PHC2z()).dot(Htemp)[2][5]))
	Htemp[:,6] *= np.sign(np.real(np.conjugate(np.transpose(Htemp)).dot(PHC2z()).dot(Htemp)[1][6]))




	HtotalEigVecs.append(np.block([[Htemp,OD2],[OD2,Htemp]]))
	HtotalEigVals.append(kin)


# REPLACE REFERS TO MY DIRECTORIES, REPLACE WITH OWN PATHS IF USING

np.savez_compressed('/n/holylfs/LABS/sachdev_lab/vectors2/HtotalEigVecsrepo_%d'%(int(sys.argv[1]),),HtotalEigVecs)
print('SHELL 1')


# Fix sign of wavefunction in different moire unit cells to enforce periodic gauge. For a model constructed with 3 shells,we find the waveufnctions are only close enough to periodic to fix the sign in the first shell of moire unit cells. However the gauge fixing condition will effect the energies only at the level of the Hartree term so we will only use the first shell where the gauge is fixed for the Hartree contribution (else the hermiticity constraints will begin to break down and the sum over form factor components which should be odd will not cancel)
for i in range(np.size(shell1momentax)):
	print('shell: ',i)

	HtotalEigVecsShell1x=[]
	for j in range(npts):
		print('pt: ',j)

		blockzer=np.full((nbands,nactive),0.0+0.0j)
		blockzerbig=np.full((nbands,nactive),0.0+0.0j)


		if (kxvals[j]+shell1momentax[i])/ktheta<0-10**(-3) and (kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta<10**(-3):

			hplus=HtotalEigVecsShell1copy[i][j][0:nbands,0:nactive]
			kxprime,kyprime=C2zrotate((kxvals[j]+shell1momentax[i]),(kyvals[j]+shell1momentay[i]))
			kxprime,kyprime=C3rotate(kxprime,kyprime)
			kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
			kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
			indexfind=np.intersect1d(kxfind,kyfind)

			if np.size(indexfind)==0:
				for k in range(np.size(shell1momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell1copy[k][indexnew][0:nbands,0:nactive]
				for k in range(np.size(shell2momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]


						U=HtotalEigVecsShell2copy[k][indexnew][0:nbands,0:nactive]
			else:
				indexnew=indexfind[0]

				U=HtotalEigVecscopy[indexnew][0:nbands,0:nactive]
			
			hminus=(C3full().dot(C3full())).dot(C2zfull()).dot(np.block([[U,blockzer],[blockzer,blockzer]]))[nbands:2*nbands,0:nactive]


			Htemp=np.block([[hplus,blockzerbig],[blockzerbig,hminus]])


		elif (kxvals[j]+shell1momentax[i])/ktheta>0-10**(-3) and -(kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta<-10**(-3):

			hplus=HtotalEigVecsShell1copy[i][j][0:nbands,0:nactive]
			kxprime,kyprime=C2zrotate((kxvals[j]+shell1momentax[i]),(kyvals[j]+shell1momentay[i]))
			kxprime,kyprime=C3rotate(kxprime,kyprime)
			kxprime,kyprime=C3rotate(kxprime,kyprime)
			kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
			kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
			indexfind=np.intersect1d(kxfind,kyfind)

			if np.size(indexfind)==0:
				for k in range(np.size(shell1momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell1copy[k][indexnew][0:nbands,0:nactive]
				for k in range(np.size(shell2momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell2copy[k][indexnew][0:nbands,0:nactive]
			else:
				indexnew=indexfind[0]

				U=HtotalEigVecscopy[indexnew][0:nbands,0:nactive]
			
			hminus=(C3full()).dot(C2zfull()).dot(np.block([[U,blockzer],[blockzer,blockzer]]))[nbands:2*nbands,0:nactive]


			Htemp=np.block([[hplus,blockzerbig],[blockzerbig,hminus]])


		elif (kxvals[j]+shell1momentax[i])/ktheta>0+10**(-3) and -(kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta>-10**(-3) and (kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta<-10**(-3):

			kxprime,kyprime=C3rotate((kxvals[j]+shell1momentax[i]),(kyvals[j]+shell1momentay[i]))

			kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
			kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
			indexfind=np.intersect1d(kxfind,kyfind)

			if np.size(indexfind)==0:
				for k in range(np.size(shell1momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell1copy[k][indexnew][0:2*nbands,0:2*nactive]
				for k in range(np.size(shell2momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell2copy[k][indexnew][0:2*nbands,0:2*nactive]
			else:
				indexnew=indexfind[0]

				U=HtotalEigVecscopy[indexnew][0:2*nbands,0:2*nactive]


			hplus=C3full().dot(C3full()).dot(U)[0:nbands,0:nactive]

			kxprime,kyprime=C2zrotate((kxvals[j]+shell1momentax[i]),(kyvals[j]+shell1momentay[i]))
			kxprime,kyprime=C3rotate(kxprime,kyprime)
			kxprime,kyprime=C3rotate(kxprime,kyprime)
			kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
			kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
			indexfind=np.intersect1d(kxfind,kyfind)

			if np.size(indexfind)==0:
				for k in range(np.size(shell1momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]
						U=HtotalEigVecsShell1copy[k][indexnew][0:nbands,0:nactive]
				for k in range(np.size(shell2momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell2copy[k][indexnew][0:nbands,0:nactive]
			else:
				indexnew=indexfind[0]

				U=HtotalEigVecscopy[indexnew][0:nbands,0:nactive]
			
			hminus=(C3full()).dot(C2zfull()).dot(np.block([[U,blockzer],[blockzer,blockzer]]))[nbands:2*nbands,0:nactive]


			Htemp=np.block([[hplus,blockzerbig],[blockzerbig,hminus]])


		elif (kxvals[j]+shell1momentax[i])/ktheta>0+10**(-3) and (kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta>-10**(-3):

			kxprime,kyprime=C3rotate((kxvals[j]+shell1momentax[i]),(kyvals[j]+shell1momentay[i]))

			kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
			kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
			indexfind=np.intersect1d(kxfind,kyfind)

			if np.size(indexfind)==0:
				for k in range(np.size(shell1momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell1copy[k][indexnew][0:2*nbands,0:2*nactive]
						hplus=np.transpose(np.conjugate(C3full())).dot(U)[0:nbands,0:nactive]
				for k in range(np.size(shell2momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell2copy[k][indexnew][0:2*nbands,0:2*nactive]
			else:
				indexnew=indexfind[0]

				U=HtotalEigVecscopy[indexnew][0:2*nbands,0:2*nactive]
				
			hplus=np.transpose(np.conjugate(C3full())).dot(U)[0:nbands,0:nactive]


			kxprime,kyprime=C2zrotate((kxvals[j]+shell1momentax[i]),(kyvals[j]+shell1momentay[i]))
			kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
			kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
			indexfind=np.intersect1d(kxfind,kyfind)

			if np.size(indexfind)==0:
				for k in range(np.size(shell1momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell1copy[k][indexnew][0:nbands,0:nactive]
				for k in range(np.size(shell2momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell2copy[k][indexnew][0:nbands,0:nactive]
			else:
				indexnew=indexfind[0]

				U=HtotalEigVecscopy[indexnew][0:nbands,0:nactive]
			
			hminus=C2zfull().dot(np.block([[U,blockzer],[blockzer,blockzer]]))[nbands:2*nbands,0:nactive]


			Htemp=np.block([[hplus,blockzerbig],[blockzerbig,hminus]])



		elif (kxvals[j]+shell1momentax[i])/ktheta<0+10**(-3) and -(kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta>10**(-3):

			kxprime,kyprime=C3rotate((kxvals[j]+shell1momentax[i]),(kyvals[j]+shell1momentay[i]))
			kxprime,kyprime=C3rotate(kxprime,kyprime)

			kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
			kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
			indexfind=np.intersect1d(kxfind,kyfind)

			if np.size(indexfind)==0:
				for k in range(np.size(shell1momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell1copy[k][indexnew][0:2*nbands,0:2*nactive]
				for k in range(np.size(shell2momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						
						indexnew=indexfind[0]

						U=HtotalEigVecsShell2copy[k][indexnew][0:2*nbands,0:2*nactive]
			else:
				indexnew=indexfind[0]

				U=HtotalEigVecscopy[indexnew][0:2*nbands,0:2*nactive]
				
			hplus=C3full().dot(U)[0:nbands,0:nactive]

			kxprime,kyprime=C2zrotate((kxvals[j]+shell1momentax[i]),(kyvals[j]+shell1momentay[i]))
			kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
			kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
			indexfind=np.intersect1d(kxfind,kyfind)

			if np.size(indexfind)==0:
				for k in range(np.size(shell1momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell1copy[k][indexnew][0:nbands,0:nactive]
				for k in range(np.size(shell2momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell2copy[k][indexnew][0:nbands,0:nactive]
			else:
				indexnew=indexfind[0]

				U=HtotalEigVecscopy[indexnew][0:nbands,0:nactive]
			
			hminus=C2zfull().dot(np.block([[U,blockzer],[blockzer,blockzer]]))[nbands:2*nbands,0:nactive]


			Htemp=np.block([[hplus,blockzerbig],[blockzerbig,hminus]])


		elif (kxvals[j]+shell1momentax[i])/ktheta<0-10**(-3) and -(kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta<10**(-3) and (kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta>10**(-3):
			kxprime,kyprime=C3rotate((kxvals[j]+shell1momentax[i]),(kyvals[j]+shell1momentay[i]))
			kxprime,kyprime=C3rotate(kxprime,kyprime)

			kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
			kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
			indexfind=np.intersect1d(kxfind,kyfind)

			if np.size(indexfind)==0:
				for k in range(np.size(shell1momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						
						indexnew=indexfind[0]

						U=HtotalEigVecsShell1copy[k][indexnew][0:2*nbands,0:2*nactive]
				for k in range(np.size(shell2momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:

						indexnew=indexfind[0]

						U=HtotalEigVecsShell2copy[k][indexnew][0:2*nbands,0:2*nactive]
			else:
				indexnew=indexfind[0]

				U=HtotalEigVecscopy[indexnew][0:2*nbands,0:2*nactive]
				
			hplus=C3full().dot(U)[0:nbands,0:nactive]

			kxprime,kyprime=C2zrotate((kxvals[j]+shell1momentax[i]),(kyvals[j]+shell1momentay[i]))
			kxprime,kyprime=C3rotate(kxprime,kyprime)
			kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
			kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
			indexfind=np.intersect1d(kxfind,kyfind)

			if np.size(indexfind)==0:
				for k in range(np.size(shell1momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell1copy[k][indexnew][0:nbands,0:nactive]
				for k in range(np.size(shell2momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell2copy[k][indexnew][0:nbands,0:nactive]
			else:
				indexnew=indexfind[0]

				U=HtotalEigVecscopy[indexnew][0:nbands,0:nactive]
			
			hminus=(C3full()).dot(C3full()).dot(C2zfull()).dot(np.block([[U,blockzer],[blockzer,blockzer]]))[nbands:2*nbands,0:nactive]


			Htemp=np.block([[hplus,blockzerbig],[blockzerbig,hminus]])

		else:
			print('something bad happened')



		Htemp[:,4] *= -np.sign(np.real(np.conjugate(np.transpose(Htemp)).dot(PHC2z()).dot(Htemp)[2][5]))
		Htemp[:,7] *= np.sign(np.real(np.conjugate(np.transpose(Htemp)).dot(PHC2z()).dot(Htemp)[1][6]))



		Htemp[:,5] *= -np.sign(np.real(np.conjugate(np.transpose(Htemp)).dot(PHC2z()).dot(Htemp)[2][5]))
		Htemp[:,6] *= np.sign(np.real(np.conjugate(np.transpose(Htemp)).dot(PHC2z()).dot(Htemp)[1][6]))



		HtotalEigVecsShell1x.append(np.block([[Htemp,OD2],[OD2,Htemp]]))
	HtotalEigVecsShell1.append(HtotalEigVecsShell1x)
# REPLACE REFERS TO MY DIRECTORIES, REPLACE WITH OWN PATHS IF USING

np.savez_compressed('/n/holylfs/LABS/sachdev_lab/vectors2/HtotalEigVecsShell1repo_%d'%(int(sys.argv[1]),),HtotalEigVecsShell1)
print('SHELL 2')

for i in range(np.size(shell2momentax)):
	print('shell: ',i)

	HtotalEigVecsShell2x=[]
	for j in range(npts):
		print('pt: ',j)

		blockzer=np.full((nbands,nactive),0.0+0.0j)
		blockzerbig=np.full((nbands,nactive),0.0+0.0j)


		if (kxvals[j]+shell2momentax[i])/ktheta<0-10**(-3) and (kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta<10**(-3):

			hplus=HtotalEigVecsShell2copy[i][j][0:nbands,0:nactive]
			kxprime,kyprime=C2zrotate((kxvals[j]+shell2momentax[i]),(kyvals[j]+shell2momentay[i]))
			kxprime,kyprime=C3rotate(kxprime,kyprime)
			kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
			kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
			indexfind=np.intersect1d(kxfind,kyfind)

			if np.size(indexfind)==0:
				U=HtotalEigVecscopy[j][0:nbands,0:nactive]

				for k in range(np.size(shell1momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell1copy[k][indexnew][0:nbands,0:nactive]
				for k in range(np.size(shell2momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]


						U=HtotalEigVecsShell2copy[k][indexnew][0:nbands,0:nactive]
			else:
				indexnew=indexfind[0]

				U=HtotalEigVecscopy[indexnew][0:nbands,0:nactive]
			
			hminus=(C3full().dot(C3full())).dot(C2zfull()).dot(np.block([[U,blockzer],[blockzer,blockzer]]))[nbands:2*nbands,0:nactive]


			Htemp=np.block([[hplus,blockzerbig],[blockzerbig,hminus]])


		elif (kxvals[j]+shell2momentax[i])/ktheta>0-10**(-3) and -(kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta<-10**(-3):

			hplus=HtotalEigVecsShell2copy[i][j][0:nbands,0:nactive]
			kxprime,kyprime=C2zrotate((kxvals[j]+shell2momentax[i]),(kyvals[j]+shell2momentay[i]))
			kxprime,kyprime=C3rotate(kxprime,kyprime)
			kxprime,kyprime=C3rotate(kxprime,kyprime)
			kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
			kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
			indexfind=np.intersect1d(kxfind,kyfind)

			if np.size(indexfind)==0:

				U=HtotalEigVecscopy[j][0:nbands,0:nactive]

				for k in range(np.size(shell1momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell1copy[k][indexnew][0:nbands,0:nactive]
				for k in range(np.size(shell2momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell2copy[k][indexnew][0:nbands,0:nactive]

			else:
				indexnew=indexfind[0]

				U=HtotalEigVecscopy[indexnew][0:nbands,0:nactive]
			
			hminus=(C3full()).dot(C2zfull()).dot(np.block([[U,blockzer],[blockzer,blockzer]]))[nbands:2*nbands,0:nactive]


			Htemp=np.block([[hplus,blockzerbig],[blockzerbig,hminus]])


		elif (kxvals[j]+shell2momentax[i])/ktheta>0+10**(-3) and -(kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta>-10**(-3) and (kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta<-10**(-3):

			kxprime,kyprime=C3rotate((kxvals[j]+shell2momentax[i]),(kyvals[j]+shell2momentay[i]))

			kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
			kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
			indexfind=np.intersect1d(kxfind,kyfind)

			if np.size(indexfind)==0:

				U=HtotalEigVecscopy[j][0:2*nbands,0:2*nactive]

				for k in range(np.size(shell1momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]
						U=HtotalEigVecsShell1copy[k][indexnew][0:2*nbands,0:2*nactive]
				for k in range(np.size(shell2momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell2copy[k][indexnew][0:2*nbands,0:2*nactive]

			else:
				indexnew=indexfind[0]

				U=HtotalEigVecscopy[indexnew][0:2*nbands,0:2*nactive]


			hplus=C3full().dot(C3full()).dot(U)[0:nbands,0:nactive]

			kxprime,kyprime=C2zrotate((kxvals[j]+shell2momentax[i]),(kyvals[j]+shell2momentay[i]))
			kxprime,kyprime=C3rotate(kxprime,kyprime)
			kxprime,kyprime=C3rotate(kxprime,kyprime)
			kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
			kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
			indexfind=np.intersect1d(kxfind,kyfind)

			if np.size(indexfind)==0:

				U=HtotalEigVecscopy[j][0:nbands,0:nactive]

				for k in range(np.size(shell1momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell1copy[k][indexnew][0:nbands,0:nactive]
				for k in range(np.size(shell2momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell2copy[k][indexnew][0:nbands,0:nactive]

			else:
				indexnew=indexfind[0]

				U=HtotalEigVecscopy[indexnew][0:nbands,0:nactive]
			
			hminus=(C3full()).dot(C2zfull()).dot(np.block([[U,blockzer],[blockzer,blockzer]]))[nbands:2*nbands,0:nactive]


			Htemp=np.block([[hplus,blockzerbig],[blockzerbig,hminus]])


		elif (kxvals[j]+shell2momentax[i])/ktheta>0+10**(-3) and (kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta>-10**(-3):

			kxprime,kyprime=C3rotate((kxvals[j]+shell2momentax[i]),(kyvals[j]+shell2momentay[i]))

			kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
			kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
			indexfind=np.intersect1d(kxfind,kyfind)

			if np.size(indexfind)==0:

				U=HtotalEigVecscopy[j][0:2*nbands,0:2*nactive]

				for k in range(np.size(shell1momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell1copy[k][indexnew][0:2*nbands,0:2*nactive]
						hplus=np.transpose(np.conjugate(C3full())).dot(U)[0:nbands,0:nactive]
				for k in range(np.size(shell2momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell2copy[k][indexnew][0:2*nbands,0:2*nactive]

			else:
				indexnew=indexfind[0]

				U=HtotalEigVecscopy[indexnew][0:2*nbands,0:2*nactive]
				
			hplus=np.transpose(np.conjugate(C3full())).dot(U)[0:nbands,0:nactive]


			kxprime,kyprime=C2zrotate((kxvals[j]+shell2momentax[i]),(kyvals[j]+shell2momentay[i]))
			kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
			kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
			indexfind=np.intersect1d(kxfind,kyfind)

			if np.size(indexfind)==0:

				U=HtotalEigVecscopy[j][0:nbands,0:nactive]

				for k in range(np.size(shell1momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell1copy[k][indexnew][0:nbands,0:nactive]
				for k in range(np.size(shell2momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell2copy[k][indexnew][0:nbands,0:nactive]

			else:
				indexnew=indexfind[0]

				U=HtotalEigVecscopy[indexnew][0:nbands,0:nactive]
			
			hminus=C2zfull().dot(np.block([[U,blockzer],[blockzer,blockzer]]))[nbands:2*nbands,0:nactive]


			Htemp=np.block([[hplus,blockzerbig],[blockzerbig,hminus]])



		elif (kxvals[j]+shell2momentax[i])/ktheta<0+10**(-3) and -(kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta>10**(-3):

			kxprime,kyprime=C3rotate((kxvals[j]+shell2momentax[i]),(kyvals[j]+shell2momentay[i]))
			kxprime,kyprime=C3rotate(kxprime,kyprime)

			kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
			kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
			indexfind=np.intersect1d(kxfind,kyfind)

			if np.size(indexfind)==0:

				U=HtotalEigVecscopy[j][0:2*nbands,0:2*nactive]

				for k in range(np.size(shell1momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell1copy[k][indexnew][0:2*nbands,0:2*nactive]
				for k in range(np.size(shell2momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						
						indexnew=indexfind[0]

						U=HtotalEigVecsShell2copy[k][indexnew][0:2*nbands,0:2*nactive]

			else:
				indexnew=indexfind[0]
				U=HtotalEigVecscopy[indexnew][0:2*nbands,0:2*nactive]
				
			hplus=C3full().dot(U)[0:nbands,0:nactive]

			kxprime,kyprime=C2zrotate((kxvals[j]+shell2momentax[i]),(kyvals[j]+shell2momentay[i]))
			kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
			kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
			indexfind=np.intersect1d(kxfind,kyfind)

			if np.size(indexfind)==0:

				U=HtotalEigVecscopy[j][0:nbands,0:nactive]

				for k in range(np.size(shell1momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell1copy[k][indexnew][0:nbands,0:nactive]
				for k in range(np.size(shell2momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell2copy[k][indexnew][0:nbands,0:nactive]

			else:
				indexnew=indexfind[0]

				U=HtotalEigVecscopy[indexnew][0:nbands,0:nactive]
			
			hminus=C2zfull().dot(np.block([[U,blockzer],[blockzer,blockzer]]))[nbands:2*nbands,0:nactive]


			Htemp=np.block([[hplus,blockzerbig],[blockzerbig,hminus]])


		elif (kxvals[j]+shell2momentax[i])/ktheta<0-10**(-3) and -(kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta<10**(-3) and (kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta>10**(-3):
			kxprime,kyprime=C3rotate((kxvals[j]+shell2momentax[i]),(kyvals[j]+shell2momentay[i]))
			kxprime,kyprime=C3rotate(kxprime,kyprime)

			kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
			kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
			indexfind=np.intersect1d(kxfind,kyfind)

			if np.size(indexfind)==0:

				U=HtotalEigVecscopy[j][0:2*nbands,0:2*nactive]

				for k in range(np.size(shell1momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						
						indexnew=indexfind[0]

						U=HtotalEigVecsShell1copy[k][indexnew][0:2*nbands,0:2*nactive]
				for k in range(np.size(shell2momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:

						indexnew=indexfind[0]

						U=HtotalEigVecsShell2copy[k][indexnew][0:2*nbands,0:2*nactive]

			else:
				indexnew=indexfind[0]

				U=HtotalEigVecscopy[indexnew][0:2*nbands,0:2*nactive]
				
			hplus=C3full().dot(U)[0:nbands,0:nactive]

			kxprime,kyprime=C2zrotate((kxvals[j]+shell2momentax[i]),(kyvals[j]+shell2momentay[i]))
			kxprime,kyprime=C3rotate(kxprime,kyprime)
			kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
			kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
			indexfind=np.intersect1d(kxfind,kyfind)

			if np.size(indexfind)==0:

				U=HtotalEigVecscopy[j][0:nbands,0:nactive]

				for k in range(np.size(shell1momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell1copy[k][indexnew][0:nbands,0:nactive]
				for k in range(np.size(shell2momentax)):
					kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
					kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
					indexfind=np.intersect1d(kxfind,kyfind)
					if np.size(indexfind)>0:
						indexnew=indexfind[0]

						U=HtotalEigVecsShell2copy[k][indexnew][0:nbands,0:nactive]

			else:
				indexnew=indexfind[0]

				U=HtotalEigVecscopy[indexnew][0:nbands,0:nactive]
			
			hminus=(C3full()).dot(C3full()).dot(C2zfull()).dot(np.block([[U,blockzer],[blockzer,blockzer]]))[nbands:2*nbands,0:nactive]


			Htemp=np.block([[hplus,blockzerbig],[blockzerbig,hminus]])

		else:
			print('something bad happened')


		Htemp[:,4] *= -np.sign(np.real(np.conjugate(np.transpose(Htemp)).dot(PHC2z()).dot(Htemp)[2][5]))
		Htemp[:,7] *= np.sign(np.real(np.conjugate(np.transpose(Htemp)).dot(PHC2z()).dot(Htemp)[1][6]))



		Htemp[:,5] *= -np.sign(np.real(np.conjugate(np.transpose(Htemp)).dot(PHC2z()).dot(Htemp)[2][5]))
		Htemp[:,6] *= np.sign(np.real(np.conjugate(np.transpose(Htemp)).dot(PHC2z()).dot(Htemp)[1][6]))

		HtotalEigVecsShell2x.append(np.block([[Htemp,OD2],[OD2,Htemp]]))
	HtotalEigVecsShell2.append(HtotalEigVecsShell2x)


# REPLACE REFERS TO MY DIRECTORIES, REPLACE WITH OWN PATHS IF USING

np.savez_compressed('/n/holylfs/LABS/sachdev_lab/vectors2/HtotalEigVecsShell2repo_%d'%(int(sys.argv[1]),),HtotalEigVecsShell2)




HtotalEigVecscopy=HtotalEigVecs
HtotalEigVecsShell1copy=HtotalEigVecsShell1
HtotalEigVecsShell2copy=HtotalEigVecsShell2
HtotalEigVecsShell4copy=HtotalEigVecsShell4

HtotalEigVecsnew=np.full((npts,4*nbands,4*nactive),0.0j)
HtotalEigVecsShell1new=np.full((np.size(shell1momentax),npts,4*nbands,4*nactive),0.0j)
HtotalEigVecsShell2new=np.full((np.size(shell2momentax),npts,4*nbands,4*nactive),0.0j)


HtotalEigVecsShell1newnew=np.full((np.size(shell1momentax),npts,4*nbands,4*nactive),0.0j)
HtotalEigVecsShell2newnew=np.full((np.size(shell2momentax),npts,4*nbands,4*nactive),0.0j)

print('SHELL 1')
for i in range(np.size(shell1momentax)):
	print('shell: ',i)
	HtotalEigVecsShell1x=[]
	for j in range(npts):

		countminus=0
		countplus=0
		countcheck=0

		for k in range(np.size(qvalslayer)):

			checkplusx=np.where(np.round(qxvalstot/ktheta,3)==np.round((qxvalstot[k]-shell1momentax[i])/ktheta,3))
			checkplusy=np.where(np.round(qyvalstot/ktheta,3)==np.round((qyvalstot[k]-shell1momentay[i])/ktheta,3))

			checkpluslayer=np.where(np.round(qvalslayer,3)==np.round(qvalslayer[k],3))
			intersectplus=np.intersect1d(checkpluslayer,np.intersect1d(checkplusx, checkplusy))


			checkminusx=np.where(np.round(qxvalstotminus/ktheta,3)==np.round((qxvalstotminus[k]-shell1momentax[i])/ktheta,3))
			checkminusy=np.where(np.round(qyvalstotminus/ktheta,3)==np.round((qyvalstotminus[k]-shell1momentay[i])/ktheta,3))

			checkminuslayer=np.where(np.round(qvalslayer,3)==np.round(qvalslayer[k],3))
			intersectminus=np.intersect1d(checkminuslayer,np.intersect1d(checkminusx, checkminusy))
			if np.size(intersectminus)>0  and countminus==0:
				if np.abs(HtotalEigVecscopy[j][nbands+2*intersectminus[0],4])<10**(-10):
					continue
				countminus=countminus+1



				phaserowv2up=HtotalEigVecscopy[j][nbands+2*intersectminus[0],nactive:2*nactive]/HtotalEigVecsShell1copy[i][j][nbands+2*k,nactive:2*nactive]*np.abs(HtotalEigVecsShell1copy[i][j][nbands+2*k,nactive:2*nactive])/np.abs(HtotalEigVecscopy[j][nbands+2*intersectminus[0],nactive:2*nactive])

				phaserowv2down=HtotalEigVecscopy[j][3*nbands+2*intersectminus[0],3*nactive:4*nactive]/HtotalEigVecsShell1copy[i][j][3*nbands+2*k,3*nactive:4*nactive]*np.abs(HtotalEigVecsShell1copy[i][j][3*nbands+2*k,3*nactive:4*nactive])/np.abs(HtotalEigVecscopy[j][3*nbands+2*intersectminus[0],3*nactive:4*nactive])



			if np.size(intersectplus)>0  and countplus==0:
				if np.abs(HtotalEigVecscopy[j][2*intersectplus[0],0])<10**(-10) :
					continue
				countplus=countplus+1



				phaserowv1up=HtotalEigVecscopy[j][2*intersectplus[0],0:nactive]/HtotalEigVecsShell1copy[i][j][2*k,0:nactive]*np.abs(HtotalEigVecsShell1copy[i][j][2*k,0:nactive])/np.abs(HtotalEigVecscopy[j][2*intersectplus[0],0:nactive])
				
				phaserowv1down=HtotalEigVecscopy[j][2*nbands+2*intersectplus[0],2*nactive:3*nactive]/HtotalEigVecsShell1copy[i][j][2*nbands+2*k,2*nactive:3*nactive]*np.abs(HtotalEigVecsShell1copy[i][j][2*nbands+2*k,2*nactive:3*nactive])/np.abs(HtotalEigVecscopy[j][2*nbands+2*intersectplus[0],2*nactive:3*nactive])

				
			if countplus>0 and countminus>0:
				

				phaserow=np.sign(np.real((np.hstack((phaserowv1up,phaserowv2up,phaserowv1down,phaserowv2down)))))

				countcheck=countcheck+1




				phase=np.diag(phaserow)
				


				HtotalEigVecsShell1newnew[i][j]=HtotalEigVecsShell1copy[i][j].dot(phase)



print('SHELL 2')
for i in range(np.size(shell2momentax)):
	print('shell: ',i)
	HtotalEigVecsShell2x=[]
	for j in range(npts):

		countminus=0
		countplus=0

		for k in range(np.size(qvalslayer)):


			checkplusx=np.where(np.round(qxvalstot/ktheta,3)==np.round((qxvalstot[k]-shell2momentax[i])/ktheta,3))
			checkplusy=np.where(np.round(qyvalstot/ktheta,3)==np.round((qyvalstot[k]-shell2momentay[i])/ktheta,3))

			checkpluslayer=np.where(np.round(qvalslayer,3)==np.round(qvalslayer[k],3))
			intersectplus=np.intersect1d(checkpluslayer,np.intersect1d(checkplusx, checkplusy))


			checkminusx=np.where(np.round(qxvalstotminus/ktheta,3)==np.round((qxvalstotminus[k]-shell2momentax[i])/ktheta,3))
			checkminusy=np.where(np.round(qyvalstotminus/ktheta,3)==np.round((qyvalstotminus[k]-shell2momentay[i])/ktheta,3))

			checkminuslayer=np.where(np.round(qvalslayer,3)==np.round(qvalslayer[k],3))
			intersectminus=np.intersect1d(checkminuslayer,np.intersect1d(checkminusx, checkminusy))
			if np.size(intersectminus)>0  and countminus==0:
				countminus=countminus+1



				phaserowv2up=HtotalEigVecscopy[j][nbands+2*intersectminus[0],nactive:2*nactive]/HtotalEigVecsShell2copy[i][j][nbands+2*k,nactive:2*nactive]*np.abs(HtotalEigVecsShell2copy[i][j][nbands+2*k,nactive:2*nactive])/np.abs(HtotalEigVecscopy[j][nbands+2*intersectminus[0],nactive:2*nactive])

				phaserowv2down=HtotalEigVecscopy[j][3*nbands+2*intersectminus[0],3*nactive:4*nactive]/HtotalEigVecsShell2copy[i][j][3*nbands+2*k,3*nactive:4*nactive]*np.abs(HtotalEigVecsShell2copy[i][j][3*nbands+2*k,3*nactive:4*nactive])/np.abs(HtotalEigVecscopy[j][3*nbands+2*intersectminus[0],3*nactive:4*nactive])



			if np.size(intersectplus)>0  and countplus==0:
				countplus=countplus+1



				phaserowv1up=HtotalEigVecscopy[j][2*intersectplus[0],0:nactive]/HtotalEigVecsShell2copy[i][j][2*k,0:nactive]*np.abs(HtotalEigVecsShell2copy[i][j][2*k,0:nactive])/np.abs(HtotalEigVecscopy[j][2*intersectplus[0],0:nactive])
				
				phaserowv1down=HtotalEigVecscopy[j][2*nbands+2*intersectplus[0],2*nactive:3*nactive]/HtotalEigVecsShell2copy[i][j][2*nbands+2*k,2*nactive:3*nactive]*np.abs(HtotalEigVecsShell2copy[i][j][2*nbands+2*k,2*nactive:3*nactive])/np.abs(HtotalEigVecscopy[j][2*nbands+2*intersectplus[0],2*nactive:3*nactive])

				
			if countplus>0 and countminus>0:

				phaserow=np.sign(np.real((np.hstack((phaserowv1up,phaserowv2up,phaserowv1down,phaserowv2down)))))




				phase=np.diag(phaserow)


				HtotalEigVecsShell2newnew[i][j]=HtotalEigVecsShell2copy[i][j]
# REPLACE REFERS TO MY DIRECTORIES, REPLACE WITH OWN PATHS IF USING
np.savez_compressed('/n/holylfs/LABS/sachdev_lab/vectors2/HtotalEigVecsShell1repo_%d'%(int(sys.argv[1]),),HtotalEigVecsShell1newnew)
np.savez_compressed('/n/holylfs/LABS/sachdev_lab/vectors2/HtotalEigVecsShell2repo_%d'%(int(sys.argv[1]),),HtotalEigVecsShell2newnew)
